I am trying to apply fix_shake to an angle with asymmetric bonds. For some reasons the angle fixes to another value other than fix value. Here is the part of the code I use to fix the angle:
units real
bond_coeff 1 1000.00 1.0
bond_coeff 2 1000.00 1.510
angle_coeff 1 100.0 109.47
fix 1 all shake 0.0001 20 0 b 1 2 a 1
The fixes for the side bonds of the angel works fine. I also expect angle #1 to be rigid angle with 109.47 degree. But once I run lammps with this code the angle fixes to 104 degree.
Any advise is appreciated, Sahand
I am trying to apply fix_shake to an angle with asymmetric bonds. For some
reasons the angle fixes to another value other than fix value. Here is the
part of the code I use to fix the angle:
units real
bond_coeff 1 1000.00 1.0
bond_coeff 2 1000.00 1.510
angle_coeff 1 100.0 109.47
fix 1 all shake 0.0001 20 0 b 1 2 a 1
The fixes for the side bonds of the angel works fine. I also expect angle
#1 to be rigid angle with 109.47 degree. But once I run lammps with this
code the angle fixes to 104 degree.
Any advise is appreciated, Sahand
please provide a *complete* and easy to run input deck that reproduces
this issue and proves that you are using LAMMPS correctly. your chances
that somebody will look into these kind of issues are otherwise close to
zero.
axel.
Thank you Axel for quick reply.
Here is my input:
units real
atom_style full
dimension 3
read_data data.test
pair_style lj/cut/coul/long 9.0 9.0
kspace_style pppm 1.0e-4
pair_coeff * 1 0.0 0.0
pair_coeff * 2 0.0000 0.0000
pair_coeff * 3 0.0000 0.0000
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.0
bond_coeff 2 1000.00 1.510
angle_coeff 1 100.0 109.47
special_bonds lj/coul 0.0 0.0 0.5
neighbor 2.0 bin
neigh_modify every 1 delay 10 check yes
reset_timestep 0
timestep 1.0
fix 1 all shake 0.0001 20 1 b 1 2 a 1
fix 2 all nvt temp 300.0 300.0 100.0
velocity all create 300 432567 dist uniform
run 1
print “all done”
here is the data [data saved in data.test file:
LAMMPS Atom File
3 atoms
2 bonds
1 angles
0 dihedrals
3 atom types
2 bond types
1 angle types
0.000000 100.00 xlo xhi
0.000000 100.00 ylo yhi
-100.000000 120.0 zlo zhi
Masses
1 15.994000
2 1.007940
3 28.085000
Atoms
1 1 1 -0.710000 0.000000 0.000000 8.056193
2 1 3 0.310000 0.000000 0.000000 6.546193
3 1 2 0.400000 0.277602 -0.901021 8.389507
Bonds
1 1 3 1
2 2 1 2
Angles
1 1 2 1 3
here is what I get in log file for SHAKE
SHAKE stats (type/ave/delta) on step 0
1 1 0
2 1.51 0
1 109.47 0
Memory usage per processor = 4.13451 Mbytes
Step Temp E_pair E_mol TotEng Press
0 300 -0.0073924573 0 0.88685056 -6.8631418
SHAKE stats (type/ave/delta) on step 1
1 1 0
2 1.51 0
1 104.311 0
1 991.73076 -0.0081856736 0 2.9479753 1.9813033
Thank you Axel for quick reply.
>>Here is my input:
thanks. i don't see any obvious problems, so this is likely a bug and thus
requires some digging into the code.
as a workaround you can try using one of the fix rigid/small variants
instead. unfortunately, those usually require a (sometimes significantly)
shorter time step.
axel.
i think you may have found one of the oldest bugs in LAMMPS. congratulations!
i can track it back to the beginning of the subversion repository in 2006. it probably was there before.
most people won’t encounter it, because the angle constraint is almost exclusively used on water molecules, where both bonds have the same length.
the bugfix is straightforward. you just need to apply the small change indicated below to fix_shake.cpp and it would work correctly. please let us know.
thanks,
axel.
diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp
index bf0c6d6…192e888 100644
— a/src/RIGID/fix_shake.cpp
+++ b/src/RIGID/fix_shake.cpp
@@ -406,8 +406,9 @@ void FixShake::init()
// compute the angle distance as a function of 2 bond distances
angle = force->angle->equilibrium_angle(i);
- rsq = 2.0*bond_distance[bond1_type]*bond_distance[bond2_type] *
- (1.0-cos(angle));
- const double b1 = bond_distance[bond1_type];
- const double b2 = bond_distance[bond2_type];
- rsq = b1b1 + b2b2 - 2.0b1b2*cos(angle);
angle_distance[i] = sqrt(rsq);
}
}
Hi Axel,
Thank you so much for the fix and sorry for very late reply. I was off the network for more than a week. I haven’t compiled LAMMPS before but I try to apply this fixes and will keep you updated.
Best regards,
Sahand
this is included in today’s 21Oct patch - thanks Axel!
Steve