Re: [ESPResSo-devel] Dihedral force calculation

[Axel Arnold]

Dear Muhammad,

please find attached two scripts that verify that the implementation is 
correct.

- compare_potential_and_forces checks that minus the derivate of a given 
potential is indeed the force central difference up to 3rd order.

- check_dihedral makes a table of the dihedral energy, which to me looks 
exactly as described in the User's Guide.

There are however different definitions of dihedrals around, not to 
mention improper dihedrals. Therefore, you might want to compare the 
formula you expect with the one in the User's guide. Also carefully  
compare the command syntax and the provided sketch on how the particles 
are located. For cutoff reasons, the particle you define the bond on is 
one of the central particles. To me it looks like you just change the 
order of particles in the part p2 bond p1 p3 p4.

@all: as this is not the first time we have this discussion, the scripts 
are now in the samples folder, in case some other potential is challenged.

Axel

On 08/27/2012 10:28 AM, Muhammad Anwar wrote:
> 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.


--
JP Dr. Axel Arnold           Tel: +49 711 685 67609
ICP, Universität Stuttgart   Email: address@hidden
Pfaffenwaldring 27
70569 Stuttgart, Germany

check_dihedral.tcl
Description: Text Data

compare_potential_and_forces.tcl
Description: Text Data