From: David van der Spoel Date: Thu, 19 Jan 2012 19:20:46 +0000 (+0100) Subject: Added two new bonded functions. X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=9885fc19ac8cf469f4351db6cc46713710d64b36;p=alexxy%2Fgromacs.git Added two new bonded functions. The first one implements a potential for linear angles without any discontinuities at 180 degrees. A paper is being written about this. The potential has [ angles ] type 9 or F_LINEAR_ANGLES internally. A free energy option has been implemented as well for the linear angles. The second is an anharmonic polarization function, based on the normal harmonic term plus a forth order term, as introduced by the Roux group in J. Chem. Theor. Comput. 6 (2010) p. 774-786. The potential has [ polarization ] type 2 or F_ANHARM_POL internally. Change-Id: I0497ed20baf1c1c0b2a6c2354eed77bc9ed75c6d --- diff --git a/include/bondf.h b/include/bondf.h index 9d40a3fa92..db28fd376f 100644 --- a/include/bondf.h +++ b/include/bondf.h @@ -138,10 +138,10 @@ void do_dih_fup(int i,int j,int k,int l,real ddphi, * *************************************************************************/ t_ifunc bonds,g96bonds,morse_bonds,cubic_bonds,FENE_bonds,restraint_bonds; - t_ifunc angles,g96angles,cross_bond_bond,cross_bond_angle,urey_bradley,quartic_angles; + t_ifunc angles,g96angles,cross_bond_bond,cross_bond_angle,urey_bradley,quartic_angles,linear_angles; t_ifunc pdihs,idihs,rbdihs; t_ifunc tab_bonds,tab_angles,tab_dihs; - t_ifunc polarize,water_pol,thole_pol,angres,angresz,unimplemented; + t_ifunc polarize,anharm_polarize,water_pol,thole_pol,angres,angresz,unimplemented; #ifdef __cplusplus } diff --git a/include/types/idef.h b/include/types/idef.h index efa74f8059..2018095645 100644 --- a/include/types/idef.h +++ b/include/types/idef.h @@ -67,6 +67,7 @@ enum { F_RESTRBONDS, F_ANGLES, F_G96ANGLES, + F_LINEAR_ANGLES, F_CROSS_BOND_BONDS, F_CROSS_BOND_ANGLES, F_UREY_BRADLEY, @@ -101,6 +102,7 @@ enum { F_POLARIZATION, F_WATER_POL, F_THOLE_POL, + F_ANHARM_POL, F_POSRES, F_DISRES, F_DISRESVIOL, @@ -147,6 +149,7 @@ typedef union */ struct {real a,b,c; } bham; struct {real rA,krA,rB,krB; } harmonic; + struct {real klinA,aA,klinB,aB; } linangle; struct {real lowA,up1A,up2A,kA,lowB,up1B,up2B,kB; } restraint; /* No free energy supported for cubic bonds, FENE, WPOL or cross terms */ struct {real b0,kb,kcub; } cubic; @@ -156,6 +159,7 @@ typedef union struct {real theta,ktheta,r13,kUB; } u_b; struct {real theta,c[5]; } qangle; struct {real alpha; } polarize; + struct {real alpha,drcut,khyp; } anharm_polarize; struct {real al_x,al_y,al_z,rOH,rHH,rOD; } wpol; struct {real a,alpha1,alpha2,rfac; } thole; struct {real c6,c12; } lj; diff --git a/include/types/nrnb.h b/include/types/nrnb.h index 14de6c9ef4..8f18e18b5d 100644 --- a/include/types/nrnb.h +++ b/include/types/nrnb.h @@ -100,13 +100,14 @@ enum eNR_CONV, eNR_SOLVEPME,eNR_NS, eNR_RESETX, eNR_SHIFTX, eNR_CGCM, eNR_FSUM, eNR_BONDS, eNR_G96BONDS, eNR_FENEBONDS, - eNR_TABBONDS, eNR_RESTRBONDS, + eNR_TABBONDS, eNR_RESTRBONDS, eNR_LINEAR_ANGLES, eNR_ANGLES, eNR_G96ANGLES, eNR_QANGLES, eNR_TABANGLES, eNR_PROPER, eNR_IMPROPER, eNR_RB, eNR_FOURDIH, eNR_TABDIHS, eNR_DISRES, eNR_ORIRES, eNR_DIHRES, eNR_POSRES, eNR_ANGRES, eNR_ANGRESZ, eNR_MORSE, eNR_CUBICBONDS, eNR_WALLS, + eNR_POLARIZE, eNR_ANHARM_POL, eNR_WPOL, eNR_THOLE, eNR_VIRIAL, eNR_UPDATE, eNR_EXTUPDATE, eNR_STOPCM, eNR_PCOUPL, eNR_EKIN, eNR_LINCS, diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c index 273804f701..f0e5f181e3 100644 --- a/src/gmxlib/bondfree.c +++ b/src/gmxlib/bondfree.c @@ -517,6 +517,63 @@ real polarize(int nbonds, return vtot; } +real anharm_polarize(int nbonds, + const t_iatom forceatoms[],const t_iparams forceparams[], + const rvec x[],rvec f[],rvec fshift[], + const t_pbc *pbc,const t_graph *g, + real lambda,real *dvdlambda, + const t_mdatoms *md,t_fcdata *fcd, + int *global_atom_index) +{ + int i,m,ki,ai,aj,type; + real dr,dr2,fbond,vbond,fij,vtot,ksh,khyp,drcut,ddr,ddr3; + rvec dx; + ivec dt; + + vtot = 0.0; + for(i=0; (ichargeA[aj])*ONE_4PI_EPS0/forceparams[type].anharm_polarize.alpha; /* 7*/ + khyp = forceparams[type].anharm_polarize.khyp; + drcut = forceparams[type].anharm_polarize.drcut; + if (debug) + fprintf(debug,"POL: local ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh); + + ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */ + dr2 = iprod(dx,dx); /* 5 */ + dr = dr2*gmx_invsqrt(dr2); /* 10 */ + + *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond); /* 19 */ + + if (dr2 == 0.0) + continue; + + if (dr > drcut) { + ddr = dr-drcut; + ddr3 = ddr*ddr*ddr; + vbond += khyp*ddr*ddr3; + fbond -= 4*khyp*ddr3; + } + fbond *= gmx_invsqrt(dr2); /* 6 */ + vtot += vbond;/* 1*/ + + if (g) { + ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt); + ki=IVEC2IS(dt); + } + for (m=0; (mpdihs.cpB = iparams->pdihs.cpA; } break; + case F_LINEAR_ANGLES: + gmx_fio_do_real(fio,iparams->linangle.klinA); + gmx_fio_do_real(fio,iparams->linangle.aA); + gmx_fio_do_real(fio,iparams->linangle.klinB); + gmx_fio_do_real(fio,iparams->linangle.aB); + break; case F_FENEBONDS: gmx_fio_do_real(fio,iparams->fene.bm); gmx_fio_do_real(fio,iparams->fene.kb); @@ -1094,6 +1102,11 @@ void do_iparams(t_fileio *fio, t_functype ftype,t_iparams *iparams, case F_POLARIZATION: gmx_fio_do_real(fio,iparams->polarize.alpha); break; + case F_ANHARM_POL: + gmx_fio_do_real(fio,iparams->anharm_polarize.alpha); + gmx_fio_do_real(fio,iparams->anharm_polarize.drcut); + gmx_fio_do_real(fio,iparams->anharm_polarize.khyp); + break; case F_WATER_POL: if (file_version < 31) gmx_fatal(FARGS,"Old tpr files with water_polarization not supported. Make a new."); diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c index 4a362578cc..ba1ce00e95 100644 --- a/src/gmxlib/txtdump.c +++ b/src/gmxlib/txtdump.c @@ -770,6 +770,11 @@ void pr_iparams(FILE *fp,t_functype ftype,t_iparams *iparams) iparams->cross_ba.r1e,iparams->cross_ba.r2e, iparams->cross_ba.r3e,iparams->cross_ba.krt); break; + case F_LINEAR_ANGLES: + fprintf(fp,"klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e\n", + iparams->linangle.klinA,iparams->linangle.aA, + iparams->linangle.klinB,iparams->linangle.aB); + break; case F_UREY_BRADLEY: fprintf(fp,"theta=%15.8e, ktheta=%15.8e, r13=%15.8e, kUB=%15.8e\n", iparams->u_b.theta,iparams->u_b.ktheta,iparams->u_b.r13,iparams->u_b.kUB); @@ -823,6 +828,12 @@ void pr_iparams(FILE *fp,t_functype ftype,t_iparams *iparams) case F_POLARIZATION: fprintf(fp,"alpha=%15.8e\n",iparams->polarize.alpha); break; + case F_ANHARM_POL: + fprintf(fp,"alpha=%15.8e drcut=%15.8e khyp=%15.8e\n", + iparams->anharm_polarize.alpha, + iparams->anharm_polarize.drcut, + iparams->anharm_polarize.khyp); + break; case F_THOLE_POL: fprintf(fp,"a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e\n", iparams->thole.a,iparams->thole.alpha1,iparams->thole.alpha2, diff --git a/src/kernel/convparm.c b/src/kernel/convparm.c index fdef7ea815..1898c24945 100644 --- a/src/kernel/convparm.c +++ b/src/kernel/convparm.c @@ -159,6 +159,12 @@ static void assign_param(t_functype ftype,t_iparams *newparam, for(i=0; i<5; i++) newparam->qangle.c[i]=old[i+1]; break; + case F_LINEAR_ANGLES: + newparam->linangle.klinA = old[0]; + newparam->linangle.aA = old[1]; + newparam->linangle.klinB = old[2]; + newparam->linangle.aB = old[3]; + break; case F_ANGLES: case F_BONDS: case F_HARMONIC: @@ -183,6 +189,11 @@ static void assign_param(t_functype ftype,t_iparams *newparam, case F_POLARIZATION: newparam->polarize.alpha = old[0]; break; + case F_ANHARM_POL: + newparam->anharm_polarize.alpha = old[0]; + newparam->anharm_polarize.drcut = old[1]; + newparam->anharm_polarize.khyp = old[2]; + break; case F_WATER_POL: newparam->wpol.al_x =old[0]; newparam->wpol.al_y =old[1]; diff --git a/src/kernel/topdirs.c b/src/kernel/topdirs.c index 14dd59fc00..fe5e9acf50 100644 --- a/src/kernel/topdirs.c +++ b/src/kernel/topdirs.c @@ -94,6 +94,8 @@ int ifunc_index(directive d,int type) return F_QUARTIC_ANGLES; case 8: return F_TABANGLES; + case 9: + return F_LINEAR_ANGLES; default: gmx_fatal(FARGS,"Invalid angle type %d",type); break; @@ -193,7 +195,15 @@ int ifunc_index(directive d,int type) } break; case d_polarization: - return F_POLARIZATION; + switch (type) { + case 1: + return F_POLARIZATION; + case 2: + return F_ANHARM_POL; + default: + gmx_fatal(FARGS,"Invalid polarization type %d",type); + } + break; case d_thole_polarization: return F_THOLE_POL; case d_water_polarization: diff --git a/src/mdlib/shellfc.c b/src/mdlib/shellfc.c index 91d11ecd79..3295cffbf0 100644 --- a/src/mdlib/shellfc.c +++ b/src/mdlib/shellfc.c @@ -200,7 +200,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog, int i,j,nmol,type,mb,mt,a_offset,cg,mol,ftype,nra; real qS,alpha; int aS,aN=0; /* Shell and nucleus */ - int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_WATER_POL }; + int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL }; #define NBT asize(bondtypes) t_iatom *ia; gmx_mtop_atomloop_block_t aloopb; @@ -304,6 +304,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog, case F_HARMONIC: case F_CUBICBONDS: case F_POLARIZATION: + case F_ANHARM_POL: if (atom[ia[1]].ptype == eptShell) { aS = ia[1]; aN = ia[2]; @@ -361,6 +362,7 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog, shell[nsi].k += ffparams->iparams[type].cubic.kb; break; case F_POLARIZATION: + case F_ANHARM_POL: if (qS != atom[aS].qB) gmx_fatal(FARGS,"polarize can not be used with qA != qB"); shell[nsi].k += sqr(qS)*ONE_4PI_EPS0/