[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ESPResSo-devel] Dihedral force calculation
From: |
Muhammad Anwar |
Subject: |
[ESPResSo-devel] Dihedral force calculation |
Date: |
Mon, 27 Aug 2012 10:28:39 +0200 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:14.0) Gecko/20120714 Thunderbird/14.0 |
Hi dear Developers,
I am using dihedral bond potential for my simulations. I found there is
a mistake in force calculation. Forces are shuffled. The force which
should be on Particle 1, it is on part 2. The force which should be on
Part 3 it is on Part 1.
You can change this
MDINLINE int calc_dihedral_force(Particle *p2, Particle *p1, Particle
*p3, Particle *p4,
Bonded_ia_parameters *iaparams, double force2[3],
double force1[2], double force3[2])
Change this like
MDINLINE int calc_dihedral_force(Particle *p2, Particle *p1, Particle
*p3, Particle *p4,
Bonded_ia_parameters *iaparams, double force1[3],
double force2[3], double force3[3])
Then change
/* store dihedral forces */
for(i=0;i<3;i++) {
force1[i] = fac*v23Xf1[i];
force2[i] = fac*(v34Xf4[i] - v12Xf1[i] - v23Xf1[i]);
force3[i] = fac*(v12Xf1[i] - v23Xf4[i] - v34Xf4[i]);
}
with
for(i=0;i<3;i++) {
force1[i] = fac*v23Xf1[i];
force2[i] = fac*(v12Xf1[i] - v23Xf4[i] - v34Xf4[i]);
force3[i] = fac*(v34Xf4[i] - v12Xf1[i] - v23Xf1[i]);
}
I have tested this and the way it is available in espresso it increase
the stiffness of chain instead of decreasing. And after changing i
compared the Rg and Re with published results.
One more efficient can be:
Compute all forces and multiply with force fac like
/* calculate force components (directions) */
for(i=0;i<3;i++) {
f1[i] = (v23Xv34[i] - cosphi*v12Xv23[i])/(l_v12Xv23);
f4[i] = (v12Xv23[i] - cosphi*v23Xv34[i])/(l_v23Xv34);
f2[i] = (-f1[i] + (scalar(v12,v23)/sqrlen(v23)) * f1[i] -
(scalar(v34,v23)/sqrlen(v23)) * f4[i]);
f3[i] = (-f4[i] - (scalar(v12,v23)/sqrlen(v23)) * f1[i] +
(scalar(v34,v23)/sqrlen(v23))* f4[i]);
}
/* store dihedral forces */
for(i=0;i<3;i++) {
force1[i] = fac * f1[i];
force2[i] = fac * f2[i];
force3[i] = fac * f3[i];
force4[i] = fac * f4[i];
}
I did not test this second method, just gone through literature and
found this in "Molecular Dynamics
Simulation Methods Revised by H. Bekker" in Chapter 5. They claim in
this method floating point operations are reduced by 40%.
Please have a look on this.
Thanks and have a nice day.
--
Muhammad Anwar
Theory of Soft Condensed Matter
162 a, Avenue De La Faiencerie, BRB 0.05
Universite of Luxembourg
1511 Luxembourg
Tel. (+352) 46 66 44 6106
- [ESPResSo-devel] Dihedral force calculation,
Muhammad Anwar <=