Error on thermodynamics

Hi,

I am trying to simulate water in a box around a non-movable iron slab. I’ve used the tether to make sure the iron doesn’t move.

However I am not sure why the thermodynamics doesn’t work.
The Log script is attached:

LAMMPS (1 Feb 2014)
units real
atom_style full

read_data Cut_Water.lammps
orthogonal box = (-13.425 -13.277 -29.829) to (15.919 15.682 31.099)
1 by 1 by 1 MPI processor grid
reading atoms …
4067 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
velocity all create 10.0 4928459 rot yes dist gaussian

neighbor 1.5 bin

neigh_modify every 1 delay 0 check yes

pair_style lj/cut/coul/cut 12.0
#bond_style harmonic
#angle_style harmonic
#dihedral_style opls

#special_bonds lj/coul 0.0 0.0 0.5

#Assume particle 1=H, 2=O, 3=Ow, 4=Fe, 5=C, 6=N,. The initial configuration file will be written in this format.

pair_coeff 1 1 0.0460 0.400
pair_coeff 1 2 0.0836 1.7753
pair_coeff 1 3 0.0620 1.3000
pair_coeff 2 2 0.1521 3.1507
pair_coeff 2 3 0.1623 3.0550
pair_coeff 2 4 0.1127 2.6754
pair_coeff 3 3 0.1669 2.960
pair_coeff 3 4 0.1180 2.580
pair_coeff 4 4 0.08345 2.200

#bond_coeff 1 450.0 0.9572 #O-H Bond(type 1 bond) parameters, Kr and r_eq- required
#bond_coeff 2 130.0 2.58 #Fe-O bond(type 2 bond) parameter, Kr and r_eq -ZDDP paper, above Table 1

#angle_coeff 1 55 104.52 #H-O-H bond angle parameter(type1 angle), K_theta and theta_eq - required

#dihedral_coeff 1 23.1 23.2 23.3 #Assuming only one type of dihedral interaction, parameters V1, V2, V3 -required

group water type 1 2
2592 atoms in group water
group metal type 3 4
1475 atoms in group metal
#fix 1 all nve
fix 2 water rigid/nph molecule iso 0.5 0.5 1.0
1 rigid bodies with 2592 atoms
fix pull metal spring tether 50.0 0.0 0.0 0.0 25.0

dump 1 all xyz 1 dump.xyz
dump_modify 1 element H, O, O, Fe
thermo_style custom press temp etotal pe ke
thermo 1

timestep 0.0000025
run 100
Memory usage per processor = 24.0472 Mbytes
Press Temp TotEng PotEng KinEng
3.801934e+09 10.214643 1052029.6 1051984.7 44.941115
-1.8177086e+11 34.306509 6.1616814e+08 6.1616799e+08 150.93751
-nan -nan -nan 0 -nan
-nan -nan -nan 0 -nan
-nan -nan -nan 0 -nan

A time 0 pressure of 10^9 indicates a bad initial

conformation or other bad initial condittions,

Which will lead to bad dyanmics.

Steve

Hi Steve

In some simulations, I observed Step = 0 ==> Pressure = -3.12E +13… But didn’t see any nan value, and everything looks fine. Is your comment valid generally for each simulation, or only when nan values appear?

Step Temp Press Volume Density TotEng KinEng PotEng
0 417.53149 -3.12E+13 526800.62 0.99999975 -98842.994 42780.019 -141623.01
500 258.53391 -653.39994 526800.62 0.99999975 -101528.5 26489.225 -128017.72
1000 259.30452 -650.91053 526800.62 0.99999975 -101448.88 26568.181 -128017.06
1500 260.05178 -648.30889 526800.62 0.99999975 -101371.21 26644.746 -128015.95
2000 260.77627 -645.76178 526800.62 0.99999975 -101295.56 26718.976 -128014.53


Thanks

Huge pressure means huge forces, typically

from atoms too close together. If the system

can time integrate out of it, the force and pressure

will drop. If it can’t the system may blow up.

Steve

Hi,

with regards to my previous post. I’ve attached the input and the lammps script.

I’m not sure why I still have such a large pressure (30000 atm), do i need to expand the box more?

or is it an issure with my NPT?

I’m trying to simulate adsorption of the amine onto the metal.

Thanks

combine1_Jmod.lammps (473 KB)

in.OPLS_AA (4.59 KB)

Job2.o702386 (2.01 KB)

fix 1 water rigid/npt molecule temp 298.0 298.0 500.0 iso 0.5 0.5 1000.0
fix 2 amine rigid/npt molecule temp 298.0 298.0 500.0 iso 0.5 0.5 1000.0

You have these 2 lines in your input script, which is bad. Both will try

to adjust the box size. You want at most one fix npt or fix npt/rigid.

The fix rigid has a discussion near the bottom about different ways to

do a hybrid simualation where some molecules are rigid and others

are not. I don’t know if that is relevant to your model, but put all the

rigid bodies in one fix. And if you have a bunch of small molecules

then you should be using fix rigid/npt/small.

Also I suggest you monitor the volume of your system in your output.

Steve