[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ESPResSo-users] Dihedral potential
From: |
Muhammad Anwar |
Subject: |
[ESPResSo-users] Dihedral potential |
Date: |
Fri, 27 May 2011 16:33:44 +0200 |
User-agent: |
Thunderbird 2.0.0.24 (X11/20101027) |
Hi,
Dear Arnold,
I have amended the dihedral potential, I also added this potential in
the asked files. But when i try to use this potential in my TCL script,
i get error that it unknown type of potential like this
unknown bonded interaction type "DIHEDRAL3COS"
while executing
"inter 3 DIHEDRAL3COS bend1 bend2 bend3"
I am attaching the code and the snippets i added in different files if
you can have a look on it and guide me about any mistake i will really
appreciate your help.
Thank you very much for your time and have a nice weekend.
Best Regards
Muhammad Anwar
// This file is part of the ESPResSo distribution (http://www.espresso.mpg.de).
// It is therefore subject to the ESPResSo license agreement which you accepted
upon receiving the distribution
// and by which you are legally bound while utilizing this file in any form or
way.
// There is NO WARRANTY, not even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.
// You should have received a copy of that license along with this program;
// if not, refer to http://www.espresso.mpg.de/license.html where its current
version can be found, or
// write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box
3148, 55021 Mainz, Germany.
// Copyright (c) 2002-2009; all rights reserved unless otherwise stated.
#ifndef DIHEDRAL3COS_H
#define DIHEDRAL3COS_H
/** \file dihedral.h Routines to calculate the dihedral energy or/and
* and force for a particle quadruple. Note that usage of dihedrals
* increases the interaction range of bonded interactions to 2 times
* the maximal bond length! \ref forces.c
*/
#define ANGLE_NOT_DEFINED -100
/// set dihedral parameters
MDINLINE int dihedral3cos_set_params(int bond_type, double bend1, double bend2,
double bend3)
{
if(bond_type < 0)
return TCL_ERROR;
make_bond_type_exist(bond_type);
bonded_ia_params[bond_type].type = BONDED_IA_DIHEDRAL3COS;
bonded_ia_params[bond_type].num = 11;
bonded_ia_params[bond_type].p.dihedral3cos.bend1 = bend1;
bonded_ia_params[bond_type].p.dihedral3cos.bend2 = bend2;
bonded_ia_params[bond_type].p.dihedral3cos.bend3 = bend3;
mpi_bcast_ia_params(bond_type, -1);
return TCL_OK;
}
/// parse parameters for the dihedral potential
MDINLINE int inter_parse_dihedral3cos(Tcl_Interp *interp, int bond_type, int
argc, char **argv)
{
double bend1, bend2, bend3;
if (argc < 4 ) {
Tcl_AppendResult(interp, "dihedral3cos needs 3 parameters: "
"<bend1> <bend2> <bend3>", (char *) NULL);
return (TCL_ERROR);
}
if ( !ARG_IS_I(1, bend1) || !ARG_IS_D(2, bend2) || !ARG_IS_D(3, bend3) ) {
Tcl_AppendResult(interp, "dihedral3cos needs 3 parameters of types DOUBLE
DOUBLE DOUBLE: "
"<bend1> <bend2> <bend3> ", (char *) NULL);
return TCL_ERROR;
}
CHECK_VALUE(dihedral3cos_set_params(bond_type, bend1, bend2, bend3), "bond
type must be nonnegative");
}
/** Calculates the dihedral angle between particle quadruple p1, p2,
p3 and p4. The dihedral angle is the angle between the planes
specified by the particle triples (p1,p2,p3) and (p2,p3,p4).
Vectors a, b and c are the bond vectors between consequtive particles.
If the a,b or b,c are parallel the dihedral angle is not defined in which
case the routine returns phi=-1. Calling functions should check for that
(Written by: Arijit Maitra) */
MDINLINE void calc_dihedral3cos_angle(Particle *p1, Particle *p2, Particle *p3,
Particle *p4,
double a[3], double b[3], double c[3],
double aXb[3], double *l_aXb, double bXc[3],
double *l_bXc,
double *cosphi, double *phi)
{
int i;
get_mi_vector(a, p2->r.p, p1->r.p);
get_mi_vector(b, p3->r.p, p2->r.p);
get_mi_vector(c, p4->r.p, p3->r.p);
/* calculate vector product a X b and b X c */
vector_product(a, b, aXb);
vector_product(b, c, bXc);
/* calculate the unit vectors */
*l_aXb = sqrt(sqrlen(aXb));
*l_bXc = sqrt(sqrlen(bXc));
/* catch case of undefined dihedral angle */
if ( *l_aXb <= TINY_LENGTH_VALUE || *l_bXc <= TINY_LENGTH_VALUE ) { *phi =
-1.0; *cosphi = 0; return;}
for (i=0;i<3;i++) {
aXb[i] /= *l_aXb;
bXc[i] /= *l_bXc;
}
*cosphi = scalar(aXb, bXc);
if ( fabs(fabs(*cosphi)-1) < TINY_SIN_VALUE ) *cosphi = dround(*cosphi);
/* Calculate dihedral angle */
*phi = acos(*cosphi);
if( scalar(aXb, c) < 0.0 ) *phi = (2.0*PI) - *phi;
}
/** calculate dihedral force between particles p1, p2 p3 and p4
Written by Arijit Maitra, adapted to new force interface by Hanjo,
more general new dihedral form by Ana.
*/
MDINLINE int calc_dihedral3cos_force(Particle *p2, Particle *p1, Particle *p3,
Particle *p4,
Bonded_ia_parameters *iaparams, double
force2[3],
double force1[2], double force3[2])
{
int i;
/* vectors for dihedral angle calculation */
double v12[3], v23[3], v34[3], v12Xv23[3], v23Xv34[3], l_v12Xv23, l_v23Xv34;
double v23Xf1[3], v23Xf4[3], v34Xf4[3], v12Xf1[3];
/* dihedral angle, cosine of the dihedral angle */
double phi, cosphi ;
/* force factors sinmphi_sinphi */
double fac, f1[3], f4[3];
/* dihedral angle */
calc_dihedral3cos_angle(p1, p2, p3, p4, v12, v23, v34, v12Xv23, &l_v12Xv23,
v23Xv34, &l_v23Xv34, &cosphi, &phi);
/* dihedral angle not defined - force zero */
if ( phi == -1.0 ) {
for(i=0;i<3;i++) { force1[i] = 0.0; force2[i] = 0.0; force3[i] = 0.0; }
return 0;
}
/* 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;
}
vector_product(v23, f1, v23Xf1);
vector_product(v23, f4, v23Xf4);
vector_product(v34, f4, v34Xf4);
vector_product(v12, f1, v12Xf1);
if(fabs(sin(phi)) < TINY_SIN_VALUE) {
fac = 0.5 * iaparams->p.dihedral3cos.bend1*cos(phi);
fac += 2 * iaparams->p.dihedral3cos.bend2*(cos(phi) - sin(phi));
fac += 6 * iaparams->p.dihedral3cos.bend3*cos(phi);
fac -= 1.5* iaparams->p.dihedral3cos.bend3*cos(phi);
fac /= cos(phi);}
else
fac = 0.5 * iaparams->p.dihedral3cos.bend1*sin(phi);
fac += 2 * iaparams->p.dihedral3cos.bend2*cos(phi);
fac += 6 * iaparams->p.dihedral3cos.bend3*sin(phi);
fac -= 1.5* iaparams->p.dihedral3cos.bend3*sin(phi);
fac /= sin(phi);
/* 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]);
}
return 0;
}
/** calculate dihedral energy between particles p1, p2 p3 and p4
Written by Arijit Maitra, adapted to new force interface by Hanjo */
MDINLINE int dihedral3cos_energy(Particle *p1, Particle *p2, Particle *p3,
Particle *p4,
Bonded_ia_parameters *iaparams, double *_energy)
{
/* vectors for dihedral calculations. */
double v12[3], v23[3], v34[3], v12Xv23[3], v23Xv34[3], l_v12Xv23, l_v23Xv34;
/* dihedral angle, cosine of the dihedral angle */
double phi, cosphi;
/* energy factors */
double fac;
calc_dihedral3cos_angle(p1, p2, p3, p4, v12, v23, v34, v12Xv23, &l_v12Xv23,
v23Xv34, &l_v23Xv34, &cosphi, &phi);
fac = 0.5 *iaparams->p.dihedral3cos.bend1 * (1 - cosphi);
fac += 0.5 * iaparams->p.dihedral3cos.bend2 * (1 - cos(2*phi));
fac += 0.5 * iaparams->p.dihedral3cos.bend3 * (1 - cos(3*phi));
*_energy = fac;
return 0;
}
#endif
%PDF-1.4
%äüöÃ
2 0 obj
<>
stream
xÝYK«ë6ÞçWx]Hª-ÁEw]îú..Ü»éßïh$Y²}îs(Ƕìiæûæ!«ÿ¾ª8+>address@hidden|Õµg¸O¹«^ÝäIÃ?yȽèÑö|¦£Á=yÝrªÝµýp
@Ù"address@hidden|Gï^£S±÷ ~W]ýdüøLÐëoù!è© ä3¨äö¬¯PÛS~Eª®µò½¿³òkÁ]¬í¼£×Ýh®×~øýùËéþ<}ÂÍ*XXûügñyB>N`ÄË.ÙLWä|faddress@hidden"ÿ¬w¶ãnMµºZ;ûÊ4ãf »hö&k
È)í¡RF!·kQÔø';Êjâc5±, ¥'ÜJ*'aÆ$äÕr ȯb)ÊÓöËü\uɱcâSô¿GÌë<\»;½÷VmIG~öìçÎ/ynDÌ kïs
;ntjò²©/ÍI<