Dear all,
Lammps beginner here. For the same system, I have done two simulations, one with real units and other with metal units. The total energies and dynamics I observed were completely different and I am not exacltly sure where the issue is happening.
I am using version: LAMMPS (23 Jun 2022 - Update 4). As new users cannot upload attachments, I am unable to upload force field and data files.
-------------input with real units-----------------------------------------------------------
units real
atom_style full
boundary p p p
read_data model_data
include real_ff
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes one 10000 page 500000
timestep 1
fix nvt_fluid all nvt temp 293.15 293.15 $(500.0*dt)
thermo 1000
thermo_style custom step temp etotal ke pe
dump nvt_prd all atom 5000 real.lammpstrj
run 3000000
-----------------------output with real units-------------------------------------------------
Step Temp TotEng KinEng PotEng
0 0 -4760.6938 0 -4760.6938
1000 236.83346 -3793.2352 7208.5122 -11001.747
2000 272.13227 -2022.8279 8282.9038 -10305.732
3000 296.5722 -580.40824 9026.7833 -9607.1916
4000 297.69012 -738.44957 9060.8094 -9799.2589
5000 288.07391 -1357.6292 8768.1203 -10125.749
6000 293.19775 -1200.0399 8924.0748 -10124.115
7000 297.26376 -1160.074 9047.8325 -10207.906
8000 289.5423 -1587.9847 8812.8139 -10400.799
-------------input with metal units-----------------------------------------------------------
units metal
atom_style full
boundary p p p
read_data model_data
include metal_ff
neighbor 1.0 bin
neigh_modify delay 0 every 1 check yes one 10000 page 500000
timestep 0.001
fix nvt_fluid all nvt temp 293.15 293.15 $(500.0*dt)
thermo 1000
thermo_style custom step temp etotal ke pe
dump nvt_prd all atom 5000 metal.lammpstrj
run 3000000
-------------------output metal units------------------------------------------------------
Step Temp TotEng KinEng PotEng
0 0 350.68146 0 350.68146
1000 175.53211 422.85059 231.6805 191.17009
2000 208.23409 517.71494 274.84304 242.87189
3000 249.6121 618.23272 329.45686 288.77587
4000 284.92765 722.68076 376.06899 346.61177
5000 303.78364 775.401 400.95654 374.44446
6000 291.4873 745.24056 384.7269 360.51366
7000 283.24621 734.37825 373.84969 360.52857
8000 296.94592 762.83852 391.9316 370.90692
9000 292.35306 751.98928 385.8696 366.11968
10000 288.69605 741.26194 381.0428 360.21913
11000 296.04808 751.2778 390.74656 360.53125
12000 292.49987 747.24318 386.06336 361.17982
13000 289.08765 734.72109 381.55966 353.16142
14000 296.38084 756.58209 391.18577 365.39633
15000 295.08777 761.07903 389.47908 371.59996
16000 287.76686 738.27408 379.81639 358.4577
17000 296.00076 753.635 390.68411 362.9509
18000 293.85837 753.87567 387.85642 366.01925
19000 291.55036 738.22555 384.81013 353.41542
20000 296.1713 755.82023 390.9092 364.91103
Please let me know if you need more details.
Did you change the force field parameters?
There is quite a bit of a difference in the energy unit between real and metal.
You can always upload files to a shared folder on Dropbox, Google Drive, One Drive, or other similar services and then provide the link. For larger files, that is the preferred option, anyway.
Dear Axel,
Thank you for the reply. Yes I have changed the ff parameters as per the units. The files are in below link. Please let me know if it does not work.
https://drive.google.com/drive/folders/1rwjwD1duAmIckMbqS5Y4Tdvn_amr5yAD?usp=sharing
metal_ff and real_ff are metal and real forcefield data
model_data has the model
metal.in and real.in are input files.
Thank you.
It looks like you used a conversion factor from eV to kcal/mol of 23.0605478, but in LAMMPS units the factor is 23.060549. see for example this function, which is used to auto-convert some potential files (EAM, Tersoff, Stillinger-Weber) between metal and real units: https://github.com/lammps/lammps/blob/6a8ca34ce89edf490bd39d16b360275ed7c6266d/src/utils.cpp#L1501
Please also note that it would be good practice to set an initial kinetic energy via the velocity command and thus speed up equilibration significantly. Also reduce the timestep to at least 0.5fs since you are not using fix shake to constrain bonds involving hydrogen atoms to get acceptable energy conservation.
Thank you very much for suggestions.
Have a nice weekend.
I have just used the LAMMPS unit conversion factor, but still there is big difference between metal and real enegies. If you have any other suggestions, please let me know.
I have uploaded force fields file for metal units with correct conversion factor (metal_ff_lammpsconversion) to the shared drive (link above)
-------------------output metal units------------------------------------------------------
Step Temp TotEng KinEng PotEng
0 293.15 737.60292 386.92146 350.68146
1000 296.50926 731.37335 391.35526 340.01808
2000 292.26312 731.21793 385.75088 345.46705
3000 292.06058 743.79518 385.48355 358.31163
4000 293.89221 757.59981 387.90108 369.69874
5000 294.92886 748.59722 389.26933 359.32789
6000 292.49044 747.10797 386.05092 361.05706
7000 294.79395 757.43119 389.09127 368.33992
8000 292.8109 731.6412 386.47389 345.16731
9000 293.94561 740.87611 387.97156 352.90455
10000 294.6199 739.58695 388.86154 350.72541
11000 292.07892 736.79523 385.50777 351.28746
12000 293.98115 745.0911 388.01847 357.07262
13000 288.6262 730.51799 380.9506 349.56739
14000 290.13034 737.95611 382.93588 355.02023
15000 291.73388 751.50362 385.05236 366.45126
16000 294.04983 740.09625 388.10912 351.98713
17000 292.34366 750.63315 385.85719 364.77596
18000 295.45902 744.58974 389.96908 354.62066
19000 294.40485 738.49735 388.5777 349.91965
20000 291.44109 741.10331 384.6659 356.43741
Of course the energies are in different units. Did you account for that?
Yes, for metal units the energies in output file are in eV and for real units they are in kcal/mol.
When I convert the energies to compare they did not match. Also the signs of energies is opposite.
Then very likely there is an error in your conversion of the force field parameters.
Please compare to the attached log files that use on-the-fly-conversion for the real_ff file for both unit sets (using immediate variables) and output the energies in both units in the same fashion (so you can compare the files side-by-side) with custom thermo column headers.
Outside of rounding differences, the two trajectories in different units agree.
log.metal (17.0 KB)
log.real (16.6 KB)
They have:
Setting up Verlet run ...
Unit style : real
Current step : 0
Time step : 0.5
Per MPI rank memory allocation (min/avg/max) = 19.08 | 19.46 | 19.84 Mbytes
Step Temp E_tot(kcal/mol) E_tot(eV) E_kin(kcal/mol) E_kin(eV) E_pot(kcal/mol) E_pot(eV) E_pair(kcal/mol) E_pair(eV) E_mol(kcal/mol) E_mol(eV)
0 200 1326.7165 57.531868 6087.4103 263.97508 -4760.6938 -206.44321 -14877.866 -645.16531 10117.172 438.72209
100 278.87791 1222.9172 53.030708 8488.2214 368.0841 -7265.3042 -315.05339 -18114.079 -785.50078 10848.775 470.44739
200 286.4821 1344.5622 58.305733 8719.6704 378.12068 -7375.1082 -319.81494 -18157.944 -787.40292 10782.835 467.58797
300 294.31685 1328.2667 57.599093 8958.1371 388.46157 -7629.8704 -330.86248 -18134.678 -786.39403 10504.808 455.53155
400 296.42071 1312.1142 56.898653 9022.1724 391.2384 -7710.0582 -334.33975 -18265.943 -792.08624 10555.885 457.74649
500 298.0832 1337.5722 58.002617 9072.7737 393.43268 -7735.2015 -335.43007 -18482.535 -801.47855 10747.334 466.04849
and:
Setting up Verlet run ...
Unit style : metal
Current step : 0
Time step : 0.0005
Per MPI rank memory allocation (min/avg/max) = 19.08 | 19.46 | 19.84 Mbytes
Step Temp E_tot(kcal/mol) E_tot(eV) E_kin(kcal/mol) E_kin(eV) E_pot(kcal/mol) E_pot(eV) E_pair(kcal/mol) E_pair(eV) E_mol(kcal/mol) E_mol(eV)
0 200 1326.7162 57.531856 6087.41 263.97507 -4760.6938 -206.44321 -14877.866 -645.16531 10117.172 438.72209
100 278.87792 1222.917 53.030696 8488.2212 368.08409 -7265.3042 -315.05339 -18114.079 -785.50078 10848.775 470.44739
200 286.48211 1344.5622 58.305734 8719.6704 378.12068 -7375.1081 -319.81494 -18157.944 -787.40292 10782.835 467.58798
300 294.31687 1328.2667 57.599091 8958.1372 388.46158 -7629.8706 -330.86249 -18134.678 -786.39403 10504.807 455.53154
400 296.42074 1312.1144 56.898664 9022.1728 391.23842 -7710.0584 -334.33976 -18265.943 -792.08623 10555.885 457.74648
500 298.0832 1337.3548 57.993188 9072.7734 393.43267 -7735.4186 -335.43948 -18482.753 -801.48798 10747.334 466.0485
Thank you very much for your time. I learnt new way to debug things.
Have a nice day.