Ca atom problem

Hello,

I am trying to simulate some clay structures(montmorillonites) with lammps which there are some neutralizing Na ions in their input structures. when I run the simulations for this structures everything is fine and I get the correct results, but when I change the ions in the input structures from Na to Ca, I get some different results for different sizes of my periodic box.

I appreciate any help in advance.

Meysam

Hello,

I am trying to simulate some clay structures(montmorillonites) with lammps
which there are some neutralizing Na ions in their input structures. when I
run the simulations for this structures everything is fine and I get the
correct results, but when I change the ions in the input structures from Na
to Ca, I get some different results for different sizes of my periodic box.
I appreciate any help in advance.

why should Ca^2+ *not* give different results than Na^+?

on top of that, you didn't say how you model both of them. if you
research the published literature, you should quickly find, that Ca^2+
cannot be well described by a pairwise additive potential.

axel.

Thanks Axel for the response,

I meant Ca2+ gives different results for different sizes of the box of Ca structures, however we know that some properties such as energy per atom or lattice constant should be independent of the size( as we can observe in the Na+ structures). Moreover, I’m running NPT simulations, and using clayff force field to prepare the initial structures.

Meysam

Thanks Axel for the response,
I meant Ca2+ gives different results for different sizes of the box of Ca
structures, however we know that some properties such as energy per atom or
lattice constant should be independent of the size( as we can observe in the
Na+ structures). Moreover, I'm running NPT simulations, and using clayff
force field to prepare the initial structures.

is your final system neutral?

energy per atom should be the same if you are in the global minimum,
but are you sure that your system found it?

axel.

yes I think so because energies are almost constant. These are last lines of the equilibration part of my simulation:

Step PotEng KinEng TotEng Temp Press Lz
1350000 -14316056 67670.404 -14248385 301.01192 249.41909 27.585449
1351000 -14316414 67362.421 -14249051 299.64195 501.34921 27.596719
1352000 -14315873 67567.89 -14248305 300.55592 129.52009 27.579912
1353000 -14316337 67623.445 -14248714 300.80304 -82.47928 27.59294
1354000 -14316349 67455.935 -14248893 300.05792 63.029155 27.564566
1355000 -14316087 67388.171 -14248699 299.75649 -477.99689 27.604098
1356000 -14316351 67608.934 -14248742 300.73849 184.38215 27.603607
1357000 -14316343 67211.855 -14249131 298.9722 -194.28911 27.592808
1358000 -14316188 67433.329 -14248755 299.95736 174.27452 27.601404
1359000 -14315926 67357.633 -14248568 299.62065 312.06742 27.611793
1360000 -14316255 67314.447 -14248940 299.42855 -2.7256962 27.590396
1361000 -14316330 66974.396 -14249356 297.91594 -124.26487 27.559174
1362000 -14316268 67702.681 -14248565 301.1555 -61.270204 27.587705
1363000 -14316294 67308.898 -14248986 299.40387 141.81141 27.605952
1364000 -14316253 67413.799 -14248840 299.87049 -223.2923 27.568701
1365000 -14316304 67313.23 -14248991 299.42314 -93.674368 27.596274
1366000 -14315937 67486.595 -14248450 300.1943 -207.57741 27.597853
1367000 -14316186 67546.412 -14248640 300.46038 267.79662 27.568329
1368000 -14316152 67402.981 -14248749 299.82237 -16.6901 27.574223
1369000 -14315687 67215.457 -14248472 298.98823 25.495231 27.541338
1370000 -14316170 67397.235 -14248773 299.79681 -180.21142 27.560961
1371000 -14316073 67493.93 -14248579 300.22693 -73.23393 27.617295
1372000 -14315790 67459.812 -14248330 300.07517 -275.45999 27.585354
1373000 -14316131 67245.12 -14248885 299.12017 -208.12057 27.595045
1374000 -14316039 67729.352 -14248309 301.27413 133.60069 27.563152
1375000 -14316807 67620.921 -14249186 300.79181 378.68519 27.57742
1376000 -14315940 67465.377 -14248474 300.09992 148.29721 27.549176
1377000 -14316507 67258.212 -14249248 299.17841 102.27844 27.578663
1378000 -14316188 67586.409 -14248601 300.63829 199.66101 27.558599
1379000 -14316136 67730.093 -14248406 301.27743 -123.37975 27.565131
1380000 -14316683 67606.981 -14249076 300.7298 328.99467 27.526996
1381000 -14316257 67214.882 -14249042 298.98567 110.92425 27.574626
1382000 -14316093 67595.821 -14248498 300.68016 -44.475389 27.58968
1383000 -14316126 67784.249 -14248342 301.51833 -214.97802 27.547852
1384000 -14316238 67250.772 -14248988 299.14531 -98.870142 27.564431
1385000 -14316181 67311.864 -14248870 299.41706 -107.31696 27.578048
1386000 -14316100 67581.354 -14248518 300.61581 74.46447 27.599431
1387000 -14316070 67580.149 -14248489 300.61045 192.3946 27.563749
1388000 -14316287 67572.305 -14248715 300.57556 384.46495 27.574534
1389000 -14315951 67493.664 -14248458 300.22575 460.52254 27.572483
1390000 -14316220 67448.034 -14248772 300.02278 93.861666 27.58015
1391000 -14316308 67209.059 -14249099 298.95977 426.09454 27.547815
1392000 -14316439 67302.784 -14249136 299.37667 238.59816 27.576932
1393000 -14316577 67490.347 -14249087 300.21099 60.48697 27.581249
1394000 -14315893 67757.617 -14248136 301.39987 356.66898 27.538148
1395000 -14316424 67354.641 -14249069 299.60734 85.233041 27.548612
1396000 -14316266 67176.276 -14249090 298.81394 12.659963 27.554573
1397000 -14316124 67544.219 -14248580 300.45063 87.664527 27.57372
1398000 -14316320 67507.574 -14248813 300.28762 -7.6382488 27.578532
1399000 -14316239 67492.593 -14248746 300.22098 -43.761996 27.569197
1400000 -14316176 67061.653 -14249114 298.30407 312.24043 27.536691

And everything is in real units, and I used this fix:

fix 1 all npt temp 300.0 300.0 100 x 1.0 1.0 1000 y 1.0 1.0 1000 z 1.0 1.0 1000 couple xy

Meysam

That only proves that you are in one minimum, not the global minimum.

so do you suggest that I heat up the system and then quench it (but I have water molecule in my system) or use different initial positions?

Meysam

OK, maybe this has nothing to do with your actuakl problem, but I remember running into similar problems before, so I feel compelled to give some (perhaps misplaced) suggestions.

How many atoms do you have in your system? Maybe the ground-state is inhomogeneous and your box scale is comparable with the wavelength of the inhomogeneity, perhaps you could try visualizing it, or identifying fragments/RDF to see if this can be the reason.

On the other hand, it isn’t easy to find global minima when you have things like competitive interactions (and I guess many other different potentials), so probably heating it and cooling it down again won’t help much.

pablo

ps: excuse me if this has nothing to do and I made you lose time