3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
48 #include "gmx_fatal.h"
53 #include "mtop_util.h"
56 void init_dihres(FILE *fplog
,gmx_mtop_t
*mtop
,t_inputrec
*ir
,t_fcdata
*fcd
)
60 fcd
->dihre_fc
= ir
->dihre_fc
;
62 count
= gmx_mtop_ftype_count(mtop
,F_DIHRES
);
65 fprintf(fplog
,"There are %d dihedral restraints\n",count
);
69 real
ta_dihres(int nfa
,const t_iatom forceatoms
[],const t_iparams ip
[],
70 const rvec x
[],rvec f
[],rvec fshift
[],
71 const t_pbc
*pbc
,const t_graph
*g
,
72 real lambda
,real
*dvdlambda
,
73 const t_mdatoms
*md
,t_fcdata
*fcd
,
77 int ai
,aj
,ak
,al
,i
,k
,type
,typep
,label
,power
,t1
,t2
,t3
;
78 real phi0
,phi
,ddphi
,ddp
,dp
,dp2
,dphi
,kfac
,sign
,d2r
,fc
;
79 rvec r_ij
,r_kj
,r_kl
,m
,n
;
85 type
= forceatoms
[i
++];
91 phi0
= ip
[type
].dihres
.phi
*d2r
;
92 dphi
= ip
[type
].dihres
.dphi
*d2r
;
93 kfac
= ip
[type
].dihres
.kfac
*fc
;
94 power
= ip
[type
].dihres
.power
;
95 label
= ip
[type
].dihres
.label
;
97 phi
= dih_angle(x
[ai
],x
[aj
],x
[ak
],x
[al
],pbc
,r_ij
,r_kj
,r_kl
,m
,n
,
102 fprintf(debug
,"dihres[%d]: %d %d %d %d : phi=%f, dphi=%f, kfac=%f, power=%d, label=%d\n",
103 k
++,ai
,aj
,ak
,al
,phi
,dphi
,kfac
,power
,label
);
105 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
106 * force changes if we just apply a normal harmonic.
107 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
108 * This means we will never have the periodicity problem, unless
109 * the dihedral is Pi away from phiO, which is very unlikely due to
113 if (fabs(dp
) > dphi
) {
114 /* dp cannot be outside (-2*pi,2*pi) */
128 vtot
+= 0.5*kfac
*ddp
*ddp
;
131 do_dih_fup(ai
,aj
,ak
,al
,ddphi
,r_ij
,r_kj
,r_kl
,m
,n
,
132 f
,fshift
,pbc
,g
,x
,t1
,t2
,t3
); /* 112 */