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. 
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.
Sincerely,
Pinky
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.
Sincerely,
Pinky
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.
