[lammps-users] regarding the use of Amber potential in lammps

Hi all,

I am trying to use AMBER Potential in LAMMPS. But, I’m getting the following error.

Has anyone used amber potential recently?

Thanks in advance.

LAMMPS (12 Dec 2018)

units real
atom_style full

pair_style lj/cut/coul/cut 8.0 10.0
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
#improper_style harmonic
#kspace_style pppm 0.0001

read_data data.c1.inpcrd
orthogonal box = (-40 -35 -30) to (47.226 44.662 38.794)
2 by 2 by 2 MPI processor grid
reading atoms …
87 atoms
scanning bonds …
3 = max bonds/atom
scanning angles …
6 = max angles/atom
scanning dihedrals …
17 = max dihedrals/atom
reading bonds …
90 bonds
reading angles …
165 angles
reading dihedrals …
312 dihedrals
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
18 = max # of 1-4 neighbors
21 = max # of special neighbors

pair_coeff * * 100.0 2.0

neighbor 2.0 bin
neigh_modify delay 5

timestep .001

thermo_style multi
thermo 200

fix 1 all nvt temp 275.0 275.0 100.0 tchain 1
fix 2 all shake 0.0001 20 10 t 5 6 m 1.0 a 31
ERROR: Invalid angle type index for fix shake (…/fix_shake.cpp:139)
Last command: fix 2 all shake 0.0001 20 10 t 5 6 m 1.0 a 31

the error has nothing to do with the amber force field, but you choosing the correct angle type index for applying fix shake.
this is something that you need to look up in the data file.

BTW: your (nearly) 3 year old LAMMPS version isn’t exactly recent. :wink:

Thank you.

I will check it out.

Hi Dr. Kohlmeyer,

I equilibrated the model, after equilibration it looks fine to me. But, when I tried to do the shear simulation using AMBER potential, the atoms which are in between the two slabs were not vibrating very much. Since they are amorphous polymers chains, they should vibrate crazily. I am confused about where the main problem lies. Is it in the (.data file) or should I check the potential?

Please find the following (.in file) and let me know where I am getting wrong if you have some space.

units real
dimension 3
boundary p p p

neighbor 2.5 bin
neigh_modify delay 5

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
pair_style lj/cut/coul/long 10.0
pair_modify mix arithmetic
kspace_style ewald 1.0e-4
read_data minimized.CHO

Region of Cellulose pillar

region up block INF INF 253 INF INF INF units box
region low block INF INF INF 169 INF INF units box
group up region up
group low region low
group boundary union up low
group mobile subtract all boundary
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 2 all shake 0.0001 20 10 t 5 6 m 1.0 a 31

timestep 0.01

Fix rigid boundary atoms

compute new2 mobile temp
fix 2 boundary setforce 0.0 0.0 0.0
unfix 3

Apply displacement control loading

velocity up set 0.0 0.00 0.50 units box
velocity low set 0.0 0.00 -0.50 units box
velocity mobile ramp vz -0.50 0.50 y 169 253 sum yes

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

compute str all stress/atom NULL
dump 2 all custom 200 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 10000

restart 20000 *.restart

Thanks for your patience reading.


I cannot comment on the physics of your system, I don’t know the material and its properties you are simulating and neither how well the force field you chose is suitable to model it.
The fact that it is not a regular crystal does not mean it should behave like a low viscosity liquid, so I don’t understand why you claim otherwise.

What stands out is that your simulation is very short with 10,000 steps and yet you are using an extremely short timestep of 0.01fs.
So your total simulation time covers only 1ps. Not much time to do anything for any kind of molecule, particularly a large polymer.

But at any rate, I am not your adviser yet the advice you are looking for is exactly the kind you should get from your adviser.
This is a mailing list about LAMMPS and not about how to do research with MD simulations.

Thank you so much for your kind response.

I tried with the timestep 1 fs. That time it gave me the following error:
ERROR on proc 0: Bond atoms 105 106 missing on proc 0 at step 5 (…/ntopo_bond_all.cpp:63)

I checked this error on the mailing list but could not find any suitable solution.


Check again. How to determine a suitable time step has been discussed many times. It should also be explained in any serviceable MD text book, and you can consult your adviser or more experienced colleagues.
Not to mention that jumping from 0.01fs to 1fs is a big step and that there can be other reasons like a bad initial geometry that may cause the bond error. Those can be addressed more effectively than lowering the time step to an extremely small value.

In addition, it should be very well possible to run your system with a 1fs time step, if you are using fix shake correctly. However, that was commented out in the example input you provided, without constrained bonds (to hydrogen atoms) you need to use a shorter timestep.

Thank you for the clarification.