[lammps-users] Regarding the data tarnsfer from AMBER to LAMMPS

pinky_09 · January 18, 2021, 8:09pm

Hi LAMMPS Developers and users,

I have made the model using AMBER, there I used GLYCAM_06j-1 force field for carbohydrates and gaff2 force field for lignin.
Then I used the amber2lammps.py script and generated the data file so that I can run simulations on LAMMPS.

LAMMPS data file for default_name

9066 atoms
9441 bonds
17398 angles
32314 dihedrals
0 impropers

16 atom types
21 bond types
47 angle types
59 dihedral types

5.672 43.356 xlo xhi
-58.52 21.286 ylo yhi
-45.557 46.571 zlo zhi

Masses

1 12.01
2 1.008
3 16.0
4 1.008
5 16.0
6 1.008
7 12.01
8 16.0
9 12.01
10 1.008
11 12.01
12 16.0
13 1.008
14 1.008
15 16.0
16 1.008

Pair Coeffs

1 0.10939999991572773 3.399669508450741
2 0.015700000021275418 2.29317330007821
3 0.17000000010548533 3.0000123432245225
4 0.01570000009846142 2.4713530426421655
5 0.2104000002486992 3.066473387458142
6 0.029999999999999995 0.35635948725613575
7 0.08600000012835884 3.399669507944831
8 0.20999999984182244 2.9599219016446874
9 0.1078000000751509 3.397709531257257
10 0.020800000056554115 2.600176998083396
11 0.09879999999945692 3.3152123100281
12 0.0725999997220856 3.1560978003356794
13 0.020800000019737253 2.4219972552023337
14 0.016100000017917802 2.6254785220391277
15 0.09300000011594518 3.242871333503238
16 0.00469999999789031 0.5379246460854297

Bond Coeffs

1 285.0 1.46
2 553.0 0.96
3 340.0 1.09
4 320.0 1.43
5 310.0 1.52
6 340.0 1.09
7 656.0 1.25
8 220.0 1.53
9 395.7 1.086
10 378.6 1.398
11 375.9 1.097
12 250.3 1.516
13 357.5 1.37
14 563.5 0.973
15 293.4 1.423
16 375.9 1.097
17 232.5 1.538
18 284.8 1.432
19 365.6 1.364
20 368.4 1.406
21 277.6 1.485

Angle Coeffs

1 60.0 110.00004714332478
2 100.0 112.00004800047613
3 70.0 108.50004650046124
4 60.0 110.00004714332478
5 55.0 109.50004692903693
6 70.0 107.50004607188556
7 45.0 111.00004757190045
8 60.0 110.00004714332478
9 45.0 113.50004864333965
10 50.0 111.60004782904585
11 45.0 109.50004692903693
12 45.0 111.00004757190045
13 60.0 110.00004714332478
14 100.0 112.00004800047613
15 80.0 126.00005400053564
16 63.0 112.36004826935493
17 50.0 109.50004692903693
18 70.0 115.00004928620316
19 63.0 111.10004761475803
20 48.7 119.88005114846936
21 68.8 120.0200516668362
22 39.0 107.5800458769885
23 65.6 120.77005170178906
24 47.3 110.47004728745955
25 87.3 119.20005108622101
26 62.4 109.780046819855
27 38.8 108.46004659790978
28 62.5 110.260047369346
29 46.9 109.56004706934303
30 85.3 107.97004621602035
31 49.0 107.26004608361897
32 84.6 110.19004739664149
33 66.1 117.96005066937873
34 64.9 111.51004761858671
35 85.2 110.6200476382248
36 47.5 109.56004706934303
37 65.2 112.07004797318065
38 46.8 109.80004705760963
39 87.2 119.90005138622399
40 50.7 108.5800463055642
41 68.9 118.38005050560584
42 48.5 119.86005148367252
43 68.4 120.69005189668611
44 86.6 120.8500520798498
45 66.1 121.11005173291322
46 66.3 112.48004797700933
47 85.6 108.9500469797992

Dihedral Coeffs

1 1.08 1 1
2 1.38 1 2
3 0.96 1 3
4 -0.27 1 1
5 0.05 1 3
6 -1.1 1 1
7 0.25 1 2
8 0.18 1 3
9 -0.1 1 1
10 0.95 1 2
11 0.55 1 3
12 0.17 1 3
13 0.16 1 3
14 0.15 1 3
15 0.27 1 3
16 0.1 1 3
17 0.6 1 2
18 0.45 1 1
19 0.4 1 2
20 0.02 1 1
21 -0.725 1 2
22 0.02 1 3
23 -1.0 1 1
24 0.1 1 2
25 -0.1 1 3
26 0.0 1 1
27 0.005 1 1
28 -0.41 1 2
29 0.01 1 3
30 -0.6 1 1
31 0.45 1 2
32 0.32 1 3
33 3.625 -1 2
34 0.0 -1 2
35 0.383333333 1 3
36 0.155555556 1 3
37 0.25 1 1
38 0.0 1 3
39 0.113 1 3
40 0.02 -1 1
41 0.0 1 2
42 1.01 1 3
43 1.61 -1 2
44 0.21 1 3
45 0.166666667 1 3
46 0.245 -1 2
47 0.12 1 3
48 0.51 1 1
49 0.18 1 3
50 0.08 1 3
51 0.11 1 1
52 0.29 -1 2
53 0.13 1 3
54 0.835 -1 2
55 0.9 -1 2
56 0.795 -1 2
57 0.337 1 3
58 10.5 -1 2
59 1.1 -1 2

Atoms

1 0 1 0.38399999999999995 21.9924254 -47.3421387 44.5169052
2 0 2 0.0 22.5303396 -48.2401943 44.8241926
3 0 3 -0.47100000000000003 22.927 -46.357 44.103
4 0 1 0.22499999999999998 23.736 -46.67 42.904
5 0 4 0.0 24.312 -47.58 43.09
6 0 1 0.282 24.703 -45.506 42.648 and so on.

To run this in LAMMPS, in the LAMMPS input file the following thing I am using and could not get the desired result. Am I choosing the correct potential? Or is there any other way to make the data file?
pair_style lj/charmm/coul/charmm 8.0 10.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic

I do not understand why I am making errors. Could you please guide me how to correct it?

Any suggestions would be appreciated.

Thank you.

Sincerely,
Pinky

akohlmey · January 18, 2021, 8:30pm

First off, please keep in mind, that the amber2lammps.py script is a contributed tool and about 15 years old (some contributed tools are even older).
The last bugfix was added in November 2013, thus there is a good chance that this tool may misread current Amber file formats.
In addition, as with all contributed tools, those are not checked for correctness for all cases but rather you can only assume that they worked sufficiently well for the contributing author. This is all included in the LAMMPS distribution as a starting point for people to adapt those for their own needs, but not as the fully supported and reliable converter.

The only way to determine whether it works correctly is to first try it with several tiny systems with only very few atom/bond/angle/dihedral types. Where you can compare in detail (and manually) if all settings and parameters can be properly converted. In short, you will have to learn what the functional forms used in the amber force field are, what the atom type assignments are and how they are mapped to LAMMPS atom types, what the individual parameters and partial charges are and whether they are correctly converted. If there are issues with the conversion, you may need to correct the python script or correct your data file.

what do you mean by “could not get the desired result”?
can your Amber input file be run as expected with Amber?
have you checked whether the functional form of the styles you are using in LAMMPS are the ones used by the Amber force field?
the input you quote is incomplete. where is the rest?

Axel.

pinky_09 · January 18, 2021, 8:51pm

Hi Dr. Kohlmeyer,
Thank you so much for your response. Please find the LAMMPS input file which I have used to run the simulations in LAMMPS. I was trying to do the shear simulation, the lower and upper part was moving, but the middle atoms are not interacting in the way they are supposed to. At room temperature, they should be in vibrating manner.

#MD simulation cellulose under shear
units real
dimension 3
boundary s s s

atom_style full
neighbor 2.5 bin
neigh_modify delay 5

pair_style lj/charmm/coul/charmm 8.0 10.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic

read_data data.relax

region lower block INF INF INF -50.00 INF INF
region upper block INF INF 13.80 INF INF INF
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary
#set group lower type 2
#set group upper type 3

Assign original velocities to atoms

compute new all temp
velocity all create 298.1 487639 temp new

Set up ensemble

fix 1 all nvt temp 298.1 298.1 100.0
fix 3 all temp/rescale 10 298.1 298.1 0.01 1.0
fix 1 all shake 0.0001 20 10 21 4 19 47 3 5

fix 1 all shake 0.0001 20 10 16 5 6 m 1.0 47 31

timestep 0.001

Output temperature, energy, pressure of the system

thermo_style custom step temp etotal pxx pyy pzz
thermo_modify flush yes
thermo 100
dump 1 all xyz 10000 dump.xyz

compute str all stress/atom NULL

dump 2 all custom 50000 dump.*.stress id type x y z c_str[1] c_str[2] c_str[3] c_str[4] c_str[5] c_str[6]

run 5000

shear

Fix rigid boundary atoms

compute new2 mobile temp
fix 2 boundary setforce 0.0 0.0 0.0
unfix 3
fix 4 mobile temp/rescale 100 298.1 298.1 0.01 1.0
fix_modify 4 temp new2

Apply displacement control loading

velocity upper set 1.0 0.00 0.00 units box
velocity lower set -1.0 0.00 0.00 units box
velocity mobile ramp vx -1.0 1.0 y 13.80 -50.00 sum yes

thermo_style custom step temp etotal pxx pyy pzz pxy pxz pyz
thermo_modify flush yes
thermo 100
dump 2 all xyz 1000 dump.cellulose

restart 20000 *.restart

run 50000

Thanks again for your patience reading.
Any suggestions would be greatly appreciated.

Sincerely,
Pinky

akohlmey · January 18, 2021, 9:45pm

You didn’t address several of my questions.

There are a few things that stand out in your input:

there doesn’t seem to be any equilibration
the simulation seems to be very short (and your system seem rather small for the purpose, too)
you are using fix temp/rescale as thermostat, which is a very bad choice for almost everything
your timestep is ridiculously small for the choice of units

Axel.