Hi dear LAMMPS users,
I’m going to use OPLS-UA to simulate ethanol (CH3-CH2-O-H)
In some literature, it mentioned that the bond length(CH3-CH2 and
CH2-O and O-H) and bond angle (CCO and COH) must be fixed if using
OPLS-UA , and only the dihedral between CCO and COH is changleable
through "dihedral_style opls".
And the fix shake is not workable because the topography of ethanol
does not fit.
So how do these researchers achieved the fix of these bond length and
angle or is there other substitute to achieve this.
Thanks for your advisement
1 Like
Hi dear LAMMPS users,
I’m going to use OPLS-UA to simulate ethanol (CH3-CH2-O-H)
In some literature, it mentioned that the bond length(CH3-CH2 and
CH2-O and O-H) and bond angle (CCO and COH) must be fixed if using
OPLS-UA , and only the dihedral between CCO and COH is changleable
through "dihedral_style opls".
And the fix shake is not workable because the topography of ethanol
does not fit.
So how do these researchers achieved the fix of these bond length and
angle or is there other substitute to achieve this.
the fix shake issues in LAMMPS are due to a preference of parallel
performance over program capability.
other codes that are not written to scale to a large number of
processors, especially serial codes, can implement SHAKE for such
coupled constraints.
also, when using MC, these issues usually do not occur.
the only workaround in LAMMPS i know about, is to use very large force
constants on the bonds and angles and correspondingly reduce the time
step to make sure that time integration does not diverge.
axel.
Thanks,axel.
I am really appreciate your help.
According to your help, I have used very large force constants on the
bonds and angles. With this method, the bond length and bond angle
become very stiff, but the displacement of ethanol molecular become
very big because of the very big bond stretching and bond bending
potential (just like explosion).
To solve this problem, I use run_style respa 2 10 bond 1 pair 2 in my
input script. I want to update the bond length and angle more
frequently so the bond potential would not be so big. But I found out
that the bond length and bond angle changes a lot than "run_style
verlet"( in these two situation(verlet and respa),the bond coeffs and
angle coeffs are the same).
Can you give me some advise, I am really appreciate your help.
Best
Shusen Guo
Thanks,axel.
I am really appreciate your help.
According to your help, I have used very large force constants on the
bonds and angles. With this method, the bond length and bond angle
become very stiff, but the displacement of ethanol molecular become
very big because of the very big bond stretching and bond bending
potential (just like explosion).
To solve this problem, I use run_style respa 2 10 bond 1 pair 2 in my
input script. I want to update the bond length and angle more
frequently so the bond potential would not be so big. But I found out
that the bond length and bond angle changes a lot than "run_style
verlet"( in these two situation(verlet and respa),the bond coeffs and
angle coeffs are the same).
Can you give me some advise, I am really appreciate your help.
sounds like you didn't reduce the time step as i recommended.
axel.
Thanks,axel
Really appreciate! According to your help, I have used two kinds of
input script to make a comparison.
1.
timestep 1(fs)
run_style respa 3 1 20 bond 1 pair 2 kspace 3(in this
situation,timestep for computation of the bond、pair and kspace is
0.05fs/0.05fs/1fs)
run 50000
2.in this situation, I want to make the timestep for bond is 0.05 same
as above, and timestep for pair and kspace is 0.25 in order to reduce
the timestep as you recommended.So the input script is as follows:
timestep 0.25
run_style respa 2 5 bond 1 pair 2
run 200000(make the total simulation time same)
But I found out that:
1. the bond length and angle change is more big in situation 2. I dont
understand why this happen because the timestep of bond is same in
these two situation.
2. Besides, the explosion in situation 2 is weakened. I think the
reason for this is the lower timestep for pair and kspace. Is that
right?
thanks for your help,axel
best
shusen Guo
Thanks,axel
Really appreciate! According to your help, I have used two kinds of
input script to make a comparison.
1.
timestep 1(fs)
run_style respa 3 1 20 bond 1 pair 2 kspace 3(in this
situation,timestep for computation of the bond、pair and kspace is
0.05fs/0.05fs/1fs)
run 50000
2.in this situation, I want to make the timestep for bond is 0.05 same
as above, and timestep for pair and kspace is 0.25 in order to reduce
the timestep as you recommended.So the input script is as follows:
timestep 0.25
run_style respa 2 5 bond 1 pair 2
run 200000(make the total simulation time same)
But I found out that:
1. the bond length and angle change is more big in situation 2. I dont
understand why this happen because the timestep of bond is same in
these two situation.
2. Besides, the explosion in situation 2 is weakened. I think the
reason for this is the lower timestep for pair and kspace. Is that
right?
there are two simple rules to follow.
1) if you cannot get it to work with verlet, it won't work with respa
2) find out how short the innermost timestep has to be with verlet and
only then determine the outer timestep with respa.
that translates into the following: reduce the timestep until you get
stable integration with verlet and good energy conservation and also
vary the order of magnitude of the force constants until you have a
combination that is keeping the structure sufficiently rigid. if you
make the bonds too stiff, your timestep will be too small without much
benefit, if they are too soft, you don't have a good approximation of
rigid bonds. only then, you can start using respa and then you need to
gradually increase the timestep ratio for non-bonded for as long as
you get stable integration and sufficient energy conservation.
your approach is the reverse of that, and that is why you have this
uncoordinated mess, that looks confusing and seems to be difficult to
get out of. if you can get it to work to your satisfaction, then you
have to look for alternatives.
axel.