# problem in relaxation of crystal

Dear All

I am trying to relax my model (Hydroxyapatite: ca5(po4)oh) in 300 k. So i use npt command. i modeled one unit cell .There are some problem during relaxation steps:

1.As some atoms may lost , I should set time step less than 0.00001 which it seems low value!
2. There is no bond between Oxygen-Oxygen but when I visualize my output it can be seen there is a bond!
3.when I reduce non-bonded cutoff less than .9 , none atoms lost but whwn i chose a value of 4 as cutoff distance problems come up!

here is my code :

pair_style lj/cut 9.4
pair_coeff 1 1 0.001229498 2.941736 #CAL-CAL
pair_coeff 1 2 1.85725E-05 2.986590077 #cal-O
pair_coeff 1 3 3.65136E-05 3.201945521 #CAL-P
pair_coeff 1 4 1.22815E-05 2.016228794 #CAL-H
pair_coeff 1 5 1.26564E-05 3.016160792 #CAL-O(H)

pair_coeff 2 2 0.002611857 3.032545038 #O-O
pair_coeff 2 3 0.005134916 3.251214178 #O-P
pair_coeff 2 4 1.79003E-05 2.047252709 #O-H
pair_coeff 2 5 0.001779873 3.062570761 # O-O(H)

pair_coeff 3 3 0.010095256 3.485650996 # P-P
pair_coeff 3 4 3.5192E-05 2.194874915 #P-H
pair_coeff 3 5 0.003499235 3.28340498 #P-O(H)

pair_coeff 4 4 1.22679E-07 1.382087851 # H-H
pair_coeff 4 5 0.010095256 3.485650996 #H-O(H)

pair_coeff 5 5 0.001212911 3.092893773 # O(H)-O(H)

bond_style harmonic
bond_coeff 1 2109.698 0.94
bond_coeff 2 1064.26 1.53

angle_style harmonic
angle_coeff 1 93.812 109.5

#--------------------zeroing stress---------------------------

fix relax all npt temp 300 300 0.01 iso 0 0 1000

dump 3 all xyz 1000 temp.xyz
timestep 0.00001

thermo_style custom step temp pe etotal
thermo 10000
run 1000000
unfix relax

also i attached dump file that interested readers could visualize with VMD.
regards.

temp.xyz (608 KB)

Dear All

I am trying to relax my model (Hydroxyapatite: ca5(po4)oh) in 300 k. So i use npt command. i modeled one unit cell .There are some problem during relaxation steps:

1.As some atoms may lost , I should set time step less than 0.00001 which it seems low value!
2. There is no bond between Oxygen-Oxygen but when I visualize my output it can be seen there is a bond!
3.when I reduce non-bonded cutoff less than .9 , none atoms lost but whwn i chose a value of 4 as cutoff distance problems come up!

here is my code :

pair_style lj/cut 9.4
pair_coeff 1 1 0.001229498 2.941736 #CAL-CAL
pair_coeff 1 2 1.85725E-05 2.986590077 #cal-O
pair_coeff 1 3 3.65136E-05 3.201945521 #CAL-P
pair_coeff 1 4 1.22815E-05 2.016228794 #CAL-H
pair_coeff 1 5 1.26564E-05 3.016160792 #CAL-O(H)

pair_coeff 2 2 0.002611857 3.032545038 #O-O
pair_coeff 2 3 0.005134916 3.251214178 #O-P
pair_coeff 2 4 1.79003E-05 2.047252709 #O-H
pair_coeff 2 5 0.001779873 3.062570761 # O-O(H)

pair_coeff 3 3 0.010095256 3.485650996 # P-P
pair_coeff 3 4 3.5192E-05 2.194874915 #P-H
pair_coeff 3 5 0.003499235 3.28340498 #P-O(H)

pair_coeff 4 4 1.22679E-07 1.382087851 # H-H
pair_coeff 4 5 0.010095256 3.485650996 #H-O(H)

pair_coeff 5 5 0.001212911 3.092893773 # O(H)-O(H)

bond_style harmonic
bond_coeff 1 2109.698 0.94
bond_coeff 2 1064.26 1.53

angle_style harmonic
angle_coeff 1 93.812 109.5

#--------------------zeroing stress---------------------------

fix relax all npt temp 300 300 0.01 iso 0 0 1000

dump 3 all xyz 1000 temp.xyz
timestep 0.00001

thermo_style custom step temp pe etotal
thermo 10000
run 1000000
unfix relax

also i attached dump file that interested readers could visualize with VMD.
regards.

Alex

temp.rar (244 KB)

Dear All

I am trying to relax my model (Hydroxyapatite: ca5(po4)oh) in 300 k. So i
use npt command. i modeled one unit cell .There are some problem during
relaxation steps:

1.As some atoms may lost , I should set time step less than 0.00001 which it
seems low value!

it is, but it is also not a good idea to immediately start an MD
of a non-equilibrated system with an NPT run. this has been
discussed on this list many, many times.

2. There is no bond between Oxygen-Oxygen but when I visualize my output it
can be seen there is a bond!

you don't say *how* you visualized it. most visualizers draw bonds
based on simple distance heuristics but not your bonding topology.
you have to communicate the bonding information to the viz program
to get that visualized right. again, this is a FAQ. please remember
that people responding here do volunteer their time, so you should
try to use it for real problems and not asking questions that have