4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
27 * GRowing Old MAkes el Chrono Sweat
29 static char *SRCID_calcpot_c
= "$Id$";
45 static void c_tabpot(real tabscale
, real VFtab
[],
48 int jindex
[], int jjnr
[],
50 real facel
, real charge
[],
51 real pot
[], real shiftvec
[])
54 const real nul
= 0.000000;
57 const real one
= 1.000000;
58 const real two
= 2.000000;
59 real r1
,r1t
,fijC
,eps
,eps2
,Y
,F
,Fp
,Geps
,Heps2
,VV
,FF
;
62 /* General and coulomb stuff */
63 int ii
,k
,n
,jnr
,ii3
,nj0
,nj1
,is3
,j3
,ggid
;
64 real fxJ
,fyJ
,fzJ
,fxO
,fyO
,fzO
;
65 real ixO
,iyO
,izO
,dxO
,dyO
,dzO
;
66 real txO
,tyO
,tzO
,vcO
,fsO
,qO
,rsqO
,rinv1O
,rinv2O
;
68 real jx
,jy
,jz
,shX
,shY
,shZ
,poti
;
70 /* Outer loop (over i particles) starts here */
71 for(n
=0; (n
<nri
); n
++) {
73 /* Unpack shift vector */
76 shY
= shiftvec
[is3
+1];
77 shZ
= shiftvec
[is3
+2];
79 /* Unpack I particle */
83 /* Charge of i particle(s) divided by 4 pi eps0 */
84 qO
= facel
*charge
[ii
];
86 /* Bounds for the innerloop */
90 /* Compute shifted i position */
92 iyO
= shY
+ pos
[ii3
+1];
93 izO
= shZ
+ pos
[ii3
+2];
96 /* Inner loop (over j-particles) starts right here */
97 for(k
=nj0
; (k
<nj1
); k
++) {
99 /* Unpack neighbourlist */
102 qj
= facel
*charge
[jnr
];
107 /* First one is for oxygen, with LJ */
111 rsqO
= dxO
*dxO
+ dyO
*dyO
+ dzO
*dzO
;
113 /* Doing fast invsqrt */
114 rinv1O
= invsqrt(rsqO
);
126 Geps
= eps
*VFtab
[nnn
+2];
127 Heps2
= eps2
*VFtab
[nnn
+3];
139 static void low_calc_pot(FILE *log
,int ftype
,t_forcerec
*fr
,
140 rvec x
[],t_mdatoms
*mdatoms
,real pot
[])
145 nlist
= &fr
->nlist_sr
[eNL_QQ
];
147 nlist
= &fr
->nlist_sr
[eNL_VDW
];
149 c_tabpot(fr
->tabscale
,fr
->VFtab
,nlist
->nri
,nlist
->iinr
,
150 nlist
->shift
,nlist
->jindex
,nlist
->jjnr
,
151 x
[0],fr
->epsfac
,mdatoms
->chargeA
,pot
,fr
->shift_vec
[0]);
153 fprintf(log
,"There were %d interactions\n",nlist
->nrj
);
156 void calc_pot(FILE *logf
,t_nsborder
*nsb
,t_commrec
*cr
,t_groups
*grps
,
157 t_parm
*parm
,t_topology
*top
,rvec x
[],t_forcerec
*fr
,
158 t_mdatoms
*mdatoms
,real pot
[])
160 static bool bFirst
=TRUE
;
168 fprintf(stderr
,"Doing single force calculation...\n");
171 snew(f
, nsb
->natoms
);
175 /* Reset long range forces if necessary */
176 if (fr
->bTwinRange
) {
177 clear_rvecs(nsb
->natoms
,fr
->flr
);
178 clear_rvecs(SHIFTS
,fr
->fshift_lr
);
180 if (parm
->ir
.epc
!= epcNO
)
181 calc_shifts(parm
->box
,box_size
,fr
->shift_vec
,FALSE
);
182 put_charge_groups_in_box(stdlog
,0,top
->blocks
[ebCGS
].nr
,FALSE
,
183 parm
->box
,box_size
,&(top
->blocks
[ebCGS
]),x
,
184 fr
->shift_vec
,fr
->cg_cm
);
185 /* mk_mshift(stdlog,graph,parm->box,x);*/
186 /* Do the actual neighbour searching and if twin range electrostatics
187 * also do the calculation of long range forces and energies.
190 ns(logf
,fr
,x
,f
,parm
->box
,grps
,&(parm
->ir
.opts
),top
,mdatoms
,cr
,
191 &nrnb
,nsb
,0,lam
,&dum
);
192 for(m
=0; (m
<DIM
); m
++)
193 box_size
[m
] = parm
->box
[m
][m
];
194 for(i
=0; (i
<mdatoms
->nr
); i
++)
197 pr_rvecs(debug
,0,"x",x
,mdatoms
->nr
);
198 pr_rvecs(debug
,0,"cgcm",fr
->cg_cm
,top
->blocks
[ebCGS
].nr
);
200 /* electrostatics from any atom to atoms without LJ */
201 low_calc_pot(logf
,F_SR
,fr
,x
,mdatoms
,pot
);
202 /* electrostatics from any atom to atoms with LJ */
203 low_calc_pot(logf
,F_LJ
,fr
,x
,mdatoms
,pot
);
206 void init_calcpot(int nfile
,t_filenm fnm
[],t_topology
*top
,
207 rvec
**x
,t_parm
*parm
,t_commrec
*cr
,
208 t_graph
**graph
,t_mdatoms
**mdatoms
,
209 t_nsborder
*nsb
,t_groups
*grps
,
210 t_forcerec
**fr
,real
**pot
,
213 real t
,t0
,lam
,lam0
,SAfac
;
215 char *traj
,*xtc_traj
;
221 tensor force_vir
,shake_vir
;
224 cr
->nprocs
= 1; cr
->pid
= 0; cr
->left
= 0; cr
->right
= 1;
225 open_log(ftp2fn(efLOG
,nfile
,fnm
),cr
);
227 init_lookup_table(stdlog
);
231 init_single(stdlog
,parm
,ftp2fn(efTPX
,nfile
,fnm
),top
,x
,&v
,mdatoms
,nsb
);
232 init_md(cr
,&(parm
->ir
),&t
,&t0
,&lam
,&lam0
,&SAfac
,
233 &nrnb
,&bTYZ
,top
,-1,NULL
,&traj
,&xtc_traj
,&fp_ene
,NULL
,
234 &mdebin
,grps
,vcm
,force_vir
,shake_vir
,*mdatoms
);
235 init_groups(stdlog
,*mdatoms
,&(parm
->ir
.opts
),grps
);
237 /* Calculate intramolecular shift vectors to make molecules whole again */
238 *graph
= mk_graph(&(top
->idef
),top
->atoms
.nr
,FALSE
);
239 mk_mshift(stdlog
,*graph
,parm
->box
,*x
);
241 /* Turn off watertype optimizations, to ease coding above. */
242 parm
->ir
.solvent_opt
= -1;
244 /* Turn off twin range if appropriate */
245 parm
->ir
.rvdw
= parm
->ir
.rcoulomb
;
246 parm
->ir
.rlist
= parm
->ir
.rcoulomb
;
247 fprintf(stderr
,"Using a coulomb cut-off of %g nm\n",parm
->ir
.rcoulomb
);
249 /* Turn off free energy computation */
250 parm
->ir
.bPert
= FALSE
;
252 /* Set vanderwaals to shift, to force tables */
253 parm
->ir
.vdwtype
= evdwSHIFT
;
254 parm
->ir
.rvdw_switch
= 0.0;
256 /* Initiate forcerecord */
258 init_forcerec(stdlog
,*fr
,&(parm
->ir
),&(top
->blocks
[ebMOLS
]),cr
,
259 &(top
->blocks
[ebCGS
]),&(top
->idef
),*mdatoms
,
260 nsb
,parm
->box
,FALSE
);
262 /* Remove periodicity */
263 for(m
=0; (m
<DIM
); m
++)
264 box_size
[m
] = parm
->box
[m
][m
];
265 if (parm
->ir
.eBox
!= ebtNONE
)
266 do_pbc_first(stdlog
,parm
,box_size
,*fr
,*graph
,*x
);
268 copy_mat(parm
->box
,box
);
270 snew(*pot
,nsb
->natoms
);