?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed typo in CMakeLists for mdrun-gpu
[gromacs/qmmm-gamess-us.git] / src / mdlib / shakef.c
blob6ca839f9388c97611fda1dfcd8599e32ac4c1f59
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROwing Monsters And Cloning Shrimps
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include "sysstuff.h"
42 #include "typedefs.h"
43 #include "smalloc.h"
44 #include "pbc.h"
45 #include "txtdump.h"
46 #include "vec.h"
47 #include "nrnb.h"
48 #include "constr.h"
49
50 typedef struct gmx_shakedata
51 {
52 rvec *rij;
53 real *M2;
54 real *tt;
55 real *dist2;
56 int nalloc;
57 /* SOR stuff */
58 real delta;
59 real omega;
60 real gamma;
61 } t_gmx_shakedata;
62
63 gmx_shakedata_t shake_init()
64 {
65 gmx_shakedata_t d;
66
67 snew(d,1);
68
69 d->nalloc = 0;
70 d->rij = NULL;
71 d->M2 = NULL;
72 d->tt = NULL;
73 d->dist2 = NULL;
74
75 /* SOR initialization */
76 d->delta = 0.1;
77 d->omega = 1.0;
78 d->gamma = 1000000;
79
80 return d;
81 }
82
83 static void pv(FILE *log,char *s,rvec x)
84 {
85 int m;
86
87 fprintf(log,"%5s:",s);
88 for(m=0; (m<DIM); m++)
89 fprintf(log," %10.3f",x[m]);
90 fprintf(log,"\n");
91 fflush(log);
92 }
93
94 void cshake(atom_id iatom[],int ncon,int *nnit,int maxnit,
95 real dist2[],real xp[],real rij[],real m2[],real omega,
96 real invmass[],real tt[],real lagr[],int *nerror)
97 {
98 /*
99 * r.c. van schaik and w.f. van gunsteren
100 * eth zuerich
101 * june 1992
102 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
103 */
104 /* default should be increased! MRS 8/4/2009 */
105 const real mytol=1e-10;
106
107 int ll,i,j,i3,j3,l3;
108 int ix,iy,iz,jx,jy,jz;
109 real toler,rpij2,rrpr,tx,ty,tz,diff,acor,im,jm;
110 real xh,yh,zh,rijx,rijy,rijz;
111 real tix,tiy,tiz;
112 real tjx,tjy,tjz;
113 int nit,error,iconv,nconv;
114
115 error=0;
116 nconv=1;
117 for (nit=0; (nit<maxnit) && (nconv != 0) && (error == 0); nit++) {
118 nconv=0;
119 for(ll=0; (ll<ncon) && (error == 0); ll++) {
120 l3 = 3*ll;
121 rijx = rij[l3+XX];
122 rijy = rij[l3+YY];
123 rijz = rij[l3+ZZ];
124 i = iatom[l3+1];
125 j = iatom[l3+2];
126 i3 = 3*i;
127 j3 = 3*j;
128 ix = i3+XX;
129 iy = i3+YY;
130 iz = i3+ZZ;
131 jx = j3+XX;
132 jy = j3+YY;
133 jz = j3+ZZ;
134
135 tx = xp[ix]-xp[jx];
136 ty = xp[iy]-xp[jy];
137 tz = xp[iz]-xp[jz];
138 rpij2 = tx*tx+ty*ty+tz*tz;
139 toler = dist2[ll];
140 diff = toler-rpij2;
141
142 /* iconv is zero when the error is smaller than a bound */
143 iconv = fabs(diff)*tt[ll];
144
145 if (iconv != 0) {
146 nconv = nconv + iconv;
147 rrpr = rijx*tx+rijy*ty+rijz*tz;
148
149 if (rrpr < toler*mytol)
150 error=ll;
151 else {
152 acor = omega*diff*m2[ll]/rrpr;
153 lagr[ll] += acor;
154 xh = rijx*acor;
155 yh = rijy*acor;
156 zh = rijz*acor;
157 im = invmass[i];
158 jm = invmass[j];
159 xp[ix] += xh*im;
160 xp[iy] += yh*im;
161 xp[iz] += zh*im;
162 xp[jx] -= xh*jm;
163 xp[jy] -= yh*jm;
164 xp[jz] -= zh*jm;
165 }
166 }
167 }
168 }
169 *nnit=nit;
170 *nerror=error;
171 }
172
173 int vec_shakef(FILE *fplog,gmx_shakedata_t shaked,
174 int natoms,real invmass[],int ncon,
175 t_iparams ip[],t_iatom *iatom,
176 real tol,rvec x[],rvec prime[],real omega,
177 bool bFEP,real lambda,real lagr[],
178 real invdt,rvec *v,
179 bool bCalcVir,tensor rmdr,int econq,
180 t_vetavars *vetavar)
181 {
182 rvec *rij;
183 real *M2,*tt,*dist2;
184 int maxnit=1000;
185 int nit=0,ll,i,j,type;
186 t_iatom *ia;
187 real L1,tol2,toler;
188 real mm=0.,tmp;
189 int error=0;
190 real g,vscale,rscale,rvscale;
191
192 if (ncon > shaked->nalloc)
193 {
194 shaked->nalloc = over_alloc_dd(ncon);
195 srenew(shaked->rij,shaked->nalloc);
196 srenew(shaked->M2,shaked->nalloc);
197 srenew(shaked->tt,shaked->nalloc);
198 srenew(shaked->dist2,shaked->nalloc);
199 }
200 rij = shaked->rij;
201 M2 = shaked->M2;
202 tt = shaked->tt;
203 dist2 = shaked->dist2;
204
205 L1=1.0-lambda;
206 tol2=2.0*tol;
207 ia=iatom;
208 for(ll=0; (ll<ncon); ll++,ia+=3) {
209 type = ia[0];
210 i=ia[1];
211 j=ia[2];
212
213 mm=2*(invmass[i]+invmass[j]);
214 rij[ll][XX]=x[i][XX]-x[j][XX];
215 rij[ll][YY]=x[i][YY]-x[j][YY];
216 rij[ll][ZZ]=x[i][ZZ]-x[j][ZZ];
217 M2[ll]=1.0/mm;
218 if (bFEP)
219 toler = sqr(L1*ip[type].constr.dA + lambda*ip[type].constr.dB);
220 else
221 toler = sqr(ip[type].constr.dA);
222 dist2[ll] = toler;
223 tt[ll] = 1.0/(toler*tol2);
224 }
225
226 switch (econq) {
227 case econqCoord:
228 cshake(iatom,ncon,&nit,maxnit,dist2,prime[0],rij[0],M2,omega,invmass,tt,lagr,&error);
229 break;
230 case econqVeloc:
231 crattle(iatom,ncon,&nit,maxnit,dist2,prime[0],rij[0],M2,omega,invmass,tt,lagr,&error,invdt,vetavar);
232 break;
233 }
234
235 if (nit >= maxnit) {
236 if (fplog) {
237 fprintf(fplog,"Shake did not converge in %d steps\n",maxnit);
238 }
239 fprintf(stderr,"Shake did not converge in %d steps\n",maxnit);
240 nit=0;
241 }
242 else if (error != 0)
243 {
244 if (fplog)
245 fprintf(fplog,"Inner product between old and new vector <= 0.0!\n"
246 "constraint #%d atoms %u and %u\n",
247 error-1,iatom[3*(error-1)+1]+1,iatom[3*(error-1)+2]+1);
248 fprintf(stderr,"Inner product between old and new vector <= 0.0!\n"
249 "constraint #%d atoms %u and %u\n",
250 error-1,iatom[3*(error-1)+1]+1,iatom[3*(error-1)+2]+1);
251 nit=0;
252 }
253
254 /* Constraint virial and correct the lagrange multipliers for the length */
255
256 ia=iatom;
257
258 for(ll=0; (ll<ncon); ll++,ia+=3)
259 {
260
261 if ((econq == econqCoord) && v!=NULL)
262 {
263 /* Correct the velocities */
264 mm = lagr[ll]*invmass[ia[1]]*invdt/vetavar->rscale;
265 for(i=0; i<DIM; i++)
266 {
267 v[ia[1]][i] += mm*rij[ll][i];
268 }
269 mm = lagr[ll]*invmass[ia[2]]*invdt/vetavar->rscale;
270 for(i=0; i<DIM; i++)
271 {
272 v[ia[2]][i] -= mm*rij[ll][i];
273 }
274 /* 16 flops */
275 }
276
277 /* constraint virial */
278 if (bCalcVir)
279 {
280 if (econq == econqCoord)
281 {
282 mm = lagr[ll]/vetavar->rvscale;
283 }
284 if (econq == econqVeloc)
285 {
286 mm = lagr[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
287 }
288 for(i=0; i<DIM; i++)
289 {
290 tmp = mm*rij[ll][i];
291 for(j=0; j<DIM; j++)
292 {
293 rmdr[i][j] -= tmp*rij[ll][j];
294 }
295 }
296 /* 21 flops */
297 }
298
299 /* Correct the lagrange multipliers for the length */
300 /* (more details would be useful here . . . )*/
301
302 type = ia[0];
303 if (bFEP)
304 {
305 toler = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
306 }
307 else
308 {
309 toler = ip[type].constr.dA;
310 lagr[ll] *= toler;
311 }
312 }
313
314 return nit;
315 }
316
317 static void check_cons(FILE *log,int nc,rvec x[],rvec prime[], rvec v[],
318 t_iparams ip[],t_iatom *iatom,
319 real invmass[], int econq)
320 {
321 t_iatom *ia;
322 int ai,aj;
323 int i;
324 real d,dp;
325 rvec dx,dv;
326
327 fprintf(log,
328 " i mi j mj before after should be\n");
329 ia=iatom;
330 for(i=0; (i<nc); i++,ia+=3) {
331 ai=ia[1];
332 aj=ia[2];
333 rvec_sub(x[ai],x[aj],dx);
334 d=norm(dx);
335
336 switch (econq) {
337 case econqCoord:
338 rvec_sub(prime[ai],prime[aj],dx);
339 dp=norm(dx);
340 fprintf(log,"%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
341 ai+1,1.0/invmass[ai],
342 aj+1,1.0/invmass[aj],d,dp,ip[ia[0]].constr.dA);
343 break;
344 case econqVeloc:
345 rvec_sub(v[ai],v[aj],dv);
346 d=iprod(dx,dv);
347 rvec_sub(prime[ai],prime[aj],dv);
348 dp=iprod(dx,dv);
349 fprintf(log,"%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
350 ai+1,1.0/invmass[ai],
351 aj+1,1.0/invmass[aj],d,dp,0.);
352 break;
353 }
354 }
355 }
356
357 bool bshakef(FILE *log,gmx_shakedata_t shaked,
358 int natoms,real invmass[],int nblocks,int sblock[],
359 t_idef *idef,t_inputrec *ir,matrix box,rvec x_s[],rvec prime[],
360 t_nrnb *nrnb,real *lagr,real lambda,real *dvdlambda,
361 real invdt,rvec *v,bool bCalcVir,tensor rmdr,bool bDumpOnError,int econq,t_vetavars *vetavar)
362 {
363 t_iatom *iatoms;
364 real *lam,dt_2,dvdl;
365 int i,n0,ncons,blen,type;
366 int tnit=0,trij=0;
367
368 #ifdef DEBUG
369 fprintf(log,"nblocks=%d, sblock[0]=%d\n",nblocks,sblock[0]);
370 #endif
371
372 ncons=idef->il[F_CONSTR].nr/3;
373
374 for(i=0; i<ncons; i++)
375 lagr[i] =0;
376
377 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
378 lam = lagr;
379 for(i=0; (i<nblocks); ) {
380 blen = (sblock[i+1]-sblock[i]);
381 blen /= 3;
382 n0 = vec_shakef(log,shaked,natoms,invmass,blen,idef->iparams,
383 iatoms,ir->shake_tol,x_s,prime,shaked->omega,
384 ir->efep!=efepNO,lambda,lam,invdt,v,bCalcVir,rmdr,econq,vetavar);
385
386 #ifdef DEBUGSHAKE
387 check_cons(log,blen,x_s,prime,v,idef->iparams,iatoms,invmass,econq);
388 #endif
389
390 if (n0 == 0) {
391 if (bDumpOnError && log)
392 {
393 {
394 check_cons(log,blen,x_s,prime,v,idef->iparams,iatoms,invmass,econq);
395 }
396 }
397 return FALSE;
398 }
399 tnit += n0*blen;
400 trij += blen;
401 iatoms += 3*blen; /* Increment pointer! */
402 lam += blen;
403 i++;
404 }
405 /* only for position part? */
406 if (econq == econqCoord) {
407 if (ir->efep != efepNO) {
408 dt_2 = 1/sqr(ir->delta_t);
409 dvdl = 0;
410 for(i=0; i<ncons; i++) {
411 type = idef->il[F_CONSTR].iatoms[3*i];
412 dvdl += lagr[i]*dt_2*
413 (idef->iparams[type].constr.dB-idef->iparams[type].constr.dA);
414 }
415 *dvdlambda += dvdl;
416 }
417 }
418 #ifdef DEBUG
419 fprintf(log,"tnit: %5d omega: %10.5f\n",tnit,omega);
420 #endif
421 if (ir->bShakeSOR) {
422 if (tnit > shaked->gamma) {
423 shaked->delta *= -0.5;
424 }
425 shaked->omega += shaked->delta;
426 shaked->gamma = tnit;
427 }
428 inc_nrnb(nrnb,eNR_SHAKE,tnit);
429 inc_nrnb(nrnb,eNR_SHAKE_RIJ,trij);
430 if (v)
431 inc_nrnb(nrnb,eNR_CONSTR_V,trij*2);
432 if (bCalcVir)
433 inc_nrnb(nrnb,eNR_CONSTR_VIR,trij);
434
435 return TRUE;
436 }
437
438 void crattle(atom_id iatom[],int ncon,int *nnit,int maxnit,
439 real dist2[],real vp[],real rij[],real m2[],real omega,
440 real invmass[],real tt[],real lagr[],int *nerror,real invdt,t_vetavars *vetavar)
441 {
442 /*
443 * r.c. van schaik and w.f. van gunsteren
444 * eth zuerich
445 * june 1992
446 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
447 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
448 * second part of rattle algorithm
449 */
450
451 const real mytol=1e-10;
452
453 int ll,i,j,i3,j3,l3,ii;
454 int ix,iy,iz,jx,jy,jz;
455 real toler,rijd,vpijd, vx,vy,vz,diff,acor,xdotd,fac,im,jm,imdt,jmdt;
456 real xh,yh,zh,rijx,rijy,rijz;
457 real tix,tiy,tiz;
458 real tjx,tjy,tjz;
459 int nit,error,iconv,nconv;
460 real veta,vscale_nhc;
461
462 veta = vetavar->veta;
463 vscale_nhc = vetavar->vscale_nhc[0]; /* for now, just use the first state */
464
465 error=0;
466 nconv=1;
467 for (nit=0; (nit<maxnit) && (nconv != 0) && (error == 0); nit++) {
468 nconv=0;
469 for(ll=0; (ll<ncon) && (error == 0); ll++) {
470 l3 = 3*ll;
471 rijx = rij[l3+XX];
472 rijy = rij[l3+YY];
473 rijz = rij[l3+ZZ];
474 i = iatom[l3+1];
475 j = iatom[l3+2];
476 i3 = 3*i;
477 j3 = 3*j;
478 ix = i3+XX;
479 iy = i3+YY;
480 iz = i3+ZZ;
481 jx = j3+XX;
482 jy = j3+YY;
483 jz = j3+ZZ;
484 vx = vp[ix]-vp[jx];
485 vy = vp[iy]-vp[jy];
486 vz = vp[iz]-vp[jz];
487
488 vpijd = vx*rijx+vy*rijy+vz*rijz;
489 toler = dist2[ll];
490 /* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
491 xdotd = vpijd*vscale_nhc + veta*toler;
492
493 /* iconv is zero when the error is smaller than a bound */
494 iconv = fabs(xdotd)*(tt[ll]/invdt);
495
496 if (iconv != 0) {
497 nconv = nconv + iconv;
498 fac = omega*2.0*m2[ll]/toler;
499 acor = -fac*xdotd;
500 lagr[ll] += acor;
501
502 xh = rijx*acor;
503 yh = rijy*acor;
504 zh = rijz*acor;
505
506 im = invmass[i]/vscale_nhc;
507 jm = invmass[j]/vscale_nhc;
508
509 vp[ix] += xh*im;
510 vp[iy] += yh*im;
511 vp[iz] += zh*im;
512 vp[jx] -= xh*jm;
513 vp[jy] -= yh*jm;
514 vp[jz] -= zh*jm;
515 }
516 }
517 }
518 *nnit=nit;
519 *nerror=error;
520 }
Atempt to write QMMM interface for Gamess-US