[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* **<=**