?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / gmxlib / update.c
blob05516683bc4ede6c14a132d6e528ca71720dee8f
1 /*
2 * $Id$
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 2.0
11 *
12 * Copyright (c) 1991-1997
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
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)
20 *
21 * Also check out our WWW page:
22 * http://rugmd0.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * Gromacs Runs On Most of All Computer Systems
28 */
29 static char *SRCID_update_c = "$Id$";
30
31 #include <stdio.h>
32 #include <math.h>
33
34 #include "assert.h"
35 #include "sysstuff.h"
36 #include "smalloc.h"
37 #include "typedefs.h"
38 #include "nrnb.h"
39 #include "led.h"
40 #include "physics.h"
41 #include "invblock.h"
42 #include "macros.h"
43 #include "vveclib.h"
44 #include "vec.h"
45 #include "main.h"
46 #include "confio.h"
47 #include "update.h"
48 #include "random.h"
49 #include "futil.h"
50 #include "mshift.h"
51 #include "tgroup.h"
52 #include "force.h"
53 #include "names.h"
54 #include "txtdump.h"
55 #include "mdrun.h"
56 #include "copyrite.h"
57 #include "edsam.h"
58 #include "callf77.h"
59
60 static void calc_g(rvec x_unc,rvec x_cons,rvec g,double mdt_2)
61 {
62 int d;
63
64 for(d=0; (d<DIM);d++)
65 g[d]=(x_cons[d]-x_unc[d])*mdt_2;
66 }
67
68 static void do_shake_corr(rvec xold,rvec x,rvec v,double dt_1)
69 {
70 int d;
71
72 for(d=0; (d<DIM);d++)
73 v[d]=((double) x[d]-(double) xold[d])*dt_1;
74 }
75
76 static void do_both(rvec xold,rvec x_unc,rvec x,rvec g,
77 rvec v,real mdt_2,real dt_1)
78 {
79 real xx,yy,zz;
80
81 xx=x[XX];
82 yy=x[YY];
83 zz=x[ZZ];
84 g[XX]=(xx-x_unc[XX])*mdt_2;
85 g[YY]=(yy-x_unc[YY])*mdt_2;
86 g[ZZ]=(zz-x_unc[ZZ])*mdt_2;
87 v[XX]=(xx-xold [XX])*dt_1;
88 v[YY]=(yy-xold [YY])*dt_1;
89 v[ZZ]=(zz-xold [ZZ])*dt_1;
90 }
91
92 static void do_update(int start,int homenr,double dt,
93 rvec lamb[],t_grp_acc gstat[],
94 rvec accel[],rvec freezefac[],
95 real invmass[],ushort ptype[],
96 ushort cFREEZE[],ushort cACC[],ushort cTC[],
97 rvec x[],rvec xprime[],rvec v[],rvec vold[],rvec f[])
98 {
99 double w_dt;
100 int gf,ga,gt;
101 real vn,vv,va,vb;
102 real uold,lg;
103 int n,d;
104
105 for (n=start; (n<start+homenr); n++) {
106 w_dt = invmass[n]*dt;
107 gf = cFREEZE[n];
108 ga = cACC[n];
109 gt = cTC[n];
110
111 for (d=0; (d<DIM); d++) {
112 vn = v[n][d];
113 lg = lamb[gt][d];
114 vold[n][d] = vn;
115
116 if ((ptype[n] != eptDummy) && (ptype[n] != eptShell) && (freezefac[gf][d] != 0)) {
117 vv = lg*(vn + f[n][d]*w_dt);
118
119 /* do not scale the mean velocities u */
120 uold = gstat[ga].uold[d];
121 va = vv + accel[ga][d]*dt;
122 vb = va + (1.0-lg)*uold;
123 v[n][d] = vb;
124 xprime[n][d] = x[n][d]+vb*dt;
125 }
126 else
127 xprime[n][d] = x[n][d];
128 }
129 }
130 }
131
132 static void do_update_lang(int start,int homenr,double dt,
133 rvec freezefac[],ushort ptype[],ushort cFREEZE[],
134 rvec x[],rvec xprime[],rvec v[],rvec vold[],
135 rvec f[],real temp, real fr, int *seed)
136 {
137 const unsigned long im = 0xffff;
138 const unsigned long ia = 1093;
139 const unsigned long ic = 18257;
140 int gf;
141 real vn,vv;
142 real rfac,invfr,rhalf,jr;
143 int n,d;
144 ulong jran;
145
146 /* (r-0.5) n times: var_n = n * var_1 = n/12
147 n=4: var_n = 1/3, so multiply with 3 */
148
149 rfac = sqrt(3.0 * 2.0*BOLTZ*temp/(fr*dt));
150 rhalf = 2.0*rfac;
151 rfac = rfac/(real)im;
152 invfr = 1.0/fr;
153
154 jran = (unsigned long)((real)im*rando(seed));
155
156 for (n=start; (n<start+homenr); n++) {
157 gf = cFREEZE[n];
158 for (d=0; (d<DIM); d++) {
159 vn = v[n][d];
160 vold[n][d] = vn;
161 if ((ptype[n]!=eptDummy) && (ptype[n]!=eptShell) && freezefac[gf][d]) {
162 jran = (jran*ia+ic) & im;
163 jr = (real)jran;
164 jran = (jran*ia+ic) & im;
165 jr += (real)jran;
166 jran = (jran*ia+ic) & im;
167 jr += (real)jran;
168 jran = (jran*ia+ic) & im;
169 jr += (real)jran;
170 vv = invfr*f[n][d] + rfac * jr - rhalf;
171 v[n][d] = vv;
172 xprime[n][d] = x[n][d]+v[n][d]*dt;
173 } else
174 xprime[n][d] = x[n][d];
175 }
176 }
177 }
178
179 static void shake_calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
180 t_commrec *cr)
181 {
182 int i,m,n;
183 matrix dvir;
184
185 clear_mat(dvir);
186 for(i=0; (i<nxf); i++) {
187 for(m=0; (m<DIM); m++)
188 for(n=0; (n<DIM); n++)
189 dvir[m][n]+=x[i][m]*f[i][n];
190 }
191
192 for(m=0; (m<DIM); m++)
193 for(n=0; (n<DIM); n++)
194 vir[m][n]-=0.5*dvir[m][n];
195 }
196
197 void dump_it_all(FILE *fp,char *title,
198 int natoms,rvec x[],rvec xp[],rvec v[],
199 rvec vold[],rvec f[])
200 {
201 #ifdef DEBUG
202 fprintf(fp,"%s\n",title);
203 pr_rvecs(fp,0,"x",x,natoms);
204 pr_rvecs(fp,0,"xp",xp,natoms);
205 pr_rvecs(fp,0,"v",v,natoms);
206 pr_rvecs(fp,0,"vold",vold,natoms);
207 pr_rvecs(fp,0,"f",f,natoms);
208 #endif
209 }
210
211 void calc_ke_part(bool bFirstStep,int start,int homenr,
212 rvec vold[],rvec v[],rvec vt[],
213 t_grpopts *opts,t_mdatoms *md,t_groups *grps,
214 t_nrnb *nrnb,real lambda,real *dvdlambda)
215 {
216 int g,d,n,ga,gt;
217 rvec v_corrt;
218 real hm,vvt,vct;
219 t_grp_tcstat *tcstat=grps->tcstat;
220 t_grp_acc *grpstat=grps->grpstat;
221 real dvdl;
222
223 /* group velocities are calculated in update_grps and
224 * accumulated in acumulate_groups.
225 * Now the partial global and groups ekin.
226 */
227 for (g=0; (g<opts->ngtc); g++)
228 clear_mat(grps->tcstat[g].ekin);
229
230 if (bFirstStep) {
231 for(n=start; (n<start+homenr); n++) {
232 copy_rvec(v[n],vold[n]);
233 }
234 for (g=0; (g<opts->ngacc); g++) {
235 for(d=0; (d<DIM); d++)
236 grps->grpstat[g].ut[d]=grps->grpstat[g].u[d];
237 }
238 }
239 else {
240 for (g=0; (g<opts->ngacc); g++) {
241 for(d=0; (d<DIM); d++)
242 grps->grpstat[g].ut[d]=0.5*(grps->grpstat[g].u[d]+
243 grps->grpstat[g].uold[d]);
244 }
245 }
246
247 dvdl = 0;
248 for(n=start; (n<start+homenr); n++) {
249 ga = md->cACC[n];
250 gt = md->cTC[n];
251 hm = 0.5*md->massT[n];
252
253 for(d=0; (d<DIM); d++) {
254 vvt = 0.5*(v[n][d]+vold[n][d]);
255 vt[n][d] = vvt;
256 vct = vvt - grpstat[ga].ut[d];
257 v_corrt[d] = vct;
258 }
259 for(d=0; (d<DIM); d++) {
260 tcstat[gt].ekin[XX][d]+=hm*v_corrt[XX]*v_corrt[d];
261 tcstat[gt].ekin[YY][d]+=hm*v_corrt[YY]*v_corrt[d];
262 tcstat[gt].ekin[ZZ][d]+=hm*v_corrt[ZZ]*v_corrt[d];
263 }
264 if (md->bPerturbed[n]) {
265 dvdl-=0.5*(md->massB[n]-md->massA[n])*iprod(v_corrt,v_corrt);
266 }
267 }
268 *dvdlambda += dvdl;
269
270 #ifdef DEBUG
271 fprintf(stdlog,"ekin: U=(%12e,%12e,%12e)\n",
272 grpstat[0].ut[XX],grpstat[0].ut[YY],grpstat[0].ut[ZZ]);
273 fprintf(stdlog,"ekin: %12e\n",trace(tcstat[0].ekin));
274 #endif
275
276 inc_nrnb(nrnb,eNR_EKIN,homenr);
277 }
278
279 typedef struct {
280 atom_id iatom[3];
281 atom_id blocknr;
282 } t_sortblock;
283
284 #ifdef DEBUG
285 static void pr_sortblock(FILE *fp,char *title,int nsb,t_sortblock sb[])
286 {
287 int i;
288
289 fprintf(fp,"%s\n",title);
290 for(i=0; (i<nsb); i++)
291 fprintf(fp,"i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
292 i,sb[i].iatom[0],sb[i].iatom[1],sb[i].iatom[2],
293 sb[i].blocknr);
294 }
295 static int pcount=0;
296 #endif
297
298 int pcomp(const void *p1, const void *p2)
299 {
300 int db;
301 atom_id min1,min2,max1,max2;
302 t_sortblock *a1=(t_sortblock *)p1;
303 t_sortblock *a2=(t_sortblock *)p2;
304 #ifdef DEBUG
305 pcount++;
306 #endif
307
308 db=a1->blocknr-a2->blocknr;
309
310 if (db != 0)
311 return db;
312
313 min1=min(a1->iatom[1],a1->iatom[2]);
314 max1=max(a1->iatom[1],a1->iatom[2]);
315 min2=min(a2->iatom[1],a2->iatom[2]);
316 max2=max(a2->iatom[1],a2->iatom[2]);
317
318 if (min1 == min2)
319 return max1-max2;
320 else
321 return min1-min2;
322 }
323
324 static int icomp(const void *p1, const void *p2)
325 {
326 atom_id *a1=(atom_id *)p1;
327 atom_id *a2=(atom_id *)p2;
328
329 return (*a1)-(*a2);
330 }
331
332 void init_update(FILE *log,t_topology *top,t_inputrec *ir,
333 t_mdatoms *md,
334 int start,int homenr,
335 int *nbl,int **sbl,
336 int *nset,int **owp,int *settle_tp)
337 {
338 t_sortblock *sb;
339 t_block *blocks=&(top->blocks[ebSBLOCKS]);
340 t_idef *idef=&(top->idef);
341 t_iatom *iatom;
342 atom_id *inv_sblock;
343 int i,j,m,bnr;
344 int ncons,bstart;
345 int settle_type;
346
347 /* Output variables, initiate them right away */
348 int nblocks=0;
349 int *sblock=NULL;
350 int nsettle=0;
351 int *owptr=NULL;
352
353 if ((ir->btc) || (ir->epc != epcNO))
354 please_cite(log,"Berendsen84a");
355
356 /* Put the oxygen atoms in the owptr array */
357 nsettle=idef->il[F_SETTLE].nr/2;
358 if (nsettle > 0) {
359 snew(owptr,nsettle);
360 settle_type=idef->il[F_SETTLE].iatoms[0];
361 *settle_tp=settle_type;
362 for (j=0; (j<idef->il[F_SETTLE].nr); j+=2) {
363 if (idef->il[F_SETTLE].iatoms[j] != settle_type)
364 fatal_error(0,"More than one settle type (%d and %d)",
365 settle_type,idef->il[F_SETTLE].iatoms[j]);
366 owptr[j/2]=idef->il[F_SETTLE].iatoms[j+1];
367 #ifdef DEBUG
368 fprintf(log,"owptr[%d]=%d\n",j/2,owptr[j/2]);
369 #endif
370 }
371 /* We used to free this memory, but ED sampling needs it later on
372 * sfree(idef->il[F_SETTLE].iatoms);
373 */
374
375 please_cite(log,"Miyamoto92a");
376 }
377
378 ncons=idef->il[F_SHAKE].nr/3;
379 if (ncons > 0) {
380 bstart=(idef->pid > 0) ? blocks->multinr[idef->pid-1] : 0;
381 nblocks=blocks->multinr[idef->pid] - bstart;
382 #ifdef DEBUGIDEF
383 fprintf(stdlog,"ncons: %d, bstart: %d, nblocks: %d\n",
384 ncons,bstart,nblocks);
385 fflush(stdlog);
386 #endif
387
388 /* Calculate block number for each atom */
389 inv_sblock=make_invblock(blocks,md->nr);
390
391 /* Store the block number in temp array and
392 * sort the constraints in order of the sblock number
393 * and the atom numbers, really sorting a segment of the array!
394 */
395 #ifdef DEBUGIDEF
396 pr_idef(stdlog,0,"Before Sort",idef);
397 #endif
398 iatom=idef->il[F_SHAKE].iatoms;
399 snew(sb,ncons);
400 for(i=0; (i<ncons); i++,iatom+=3) {
401 for(m=0; (m<3); m++)
402 sb[i].iatom[m]=iatom[m];
403 sb[i].blocknr=inv_sblock[iatom[1]];
404 }
405
406 /* Now sort the blocks */
407 #ifdef DEBUG
408 pr_sortblock(log,"Before sorting",ncons,sb);
409 fprintf(log,"Going to sort constraints\n");
410 #endif
411
412 qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
413
414 #ifdef DEBUG
415 fprintf(log,"I used %d calls to pcomp\n",pcount);
416 pr_sortblock(log,"After sorting",ncons,sb);
417 #endif
418
419 iatom=idef->il[F_SHAKE].iatoms;
420 for(i=0; (i<ncons); i++,iatom+=3)
421 for(m=0; (m<DIM); m++)
422 iatom[m]=sb[i].iatom[m];
423 #ifdef DEBUGIDEF
424 pr_idef(stdlog,0,"After Sort",idef);
425 #endif
426
427 j=0;
428 snew(sblock,nblocks+1);
429 bnr=-2;
430 for(i=0; (i<ncons); i++) {
431 if (sb[i].blocknr != bnr) {
432 bnr=sb[i].blocknr;
433 sblock[j++]=3*i;
434 }
435 }
436 /* Last block... */
437 sblock[j++]=3*ncons;
438
439 if (j != (nblocks+1)) {
440 fprintf(stdlog,"bstart: %d\n",bstart);
441 fprintf(log,"j: %d, nblocks: %d, ncons: %d\n",
442 j,nblocks,ncons);
443 for(i=0; (i<ncons); i++)
444 fprintf(log,"i: %5d sb[i].blocknr: %5u\n",i,sb[i].blocknr);
445 for(j=0; (j<=nblocks); j++)
446 fprintf(log,"sblock[%3d]=%5d\n",j,(int) sblock[j]);
447 fatal_error(0,"DEATH HORROR: "
448 "top->blocks[ebSBLOCKS] does not match idef->il[F_SHAKE]");
449 }
450 sfree(sb);
451 sfree(inv_sblock);
452 }
453
454 /* Copy pointers */
455 *nbl=nblocks;
456 *sbl=sblock;
457 *nset=nsettle;
458 *owp=owptr;
459 }
460
461
462 static void init_lincs(FILE *log,t_topology *top,t_inputrec *ir,
463 t_mdatoms *md,int start,int homenr,
464 int *nbl,int **sbl,
465 int *nset,int **owp,int *settle_tp,
466 int *ncm,int *cmax,
467 rvec **r,int **bla1,int **bla2,int **blnr,int **blbnb,
468 real **bllen,real **blc,real **blcc,real **blm,
469 real **tmp1,real **tmp2,real **tmp3,
470 real **lincslam,real **bllen0,real **ddist)
471 {
472 t_idef *idef=&(top->idef);
473 t_iatom *iatom;
474 int i,j,k,n,b1,b,cen;
475 int ncons;
476 int type,a1,a2,b2,nr,n1,n2,nc4;
477 real len=0,len1,sign;
478 real im1,im2,imcen;
479
480 ncons=idef->il[F_SHAKE].nr/3;
481
482 if (ncons > 0) {
483
484 /* Make constraint-neighbour list */
485
486 snew(*r,ncons);
487 snew(*bla1,ncons);
488 snew(*bla2,ncons);
489 snew(*blnr,ncons);
490 snew(*bllen,ncons);
491 snew(*blc,ncons);
492 snew(*tmp1,ncons);
493 snew(*tmp2,ncons);
494 snew(*tmp3,ncons);
495 snew(*lincslam,ncons);
496 snew(*bllen0,ncons);
497 snew(*ddist,ncons);
498
499 iatom=idef->il[F_SHAKE].iatoms;
500
501 /* Number of coupling of a bond is defined as the number of
502 bonds directly connected to that bond (not to an atom!).
503 The constraint are divided into two groups, the first
504 group consists of bonds with 4 or less couplings, the
505 second group consists of bonds with more than 4 couplings
506 (in proteins most of the bonds have 2 or 4 couplings).
507
508 cmax: maximum number of bonds coupled to one bond
509 ncm: number of bonds with mor than 4 couplings
510 */
511
512 *cmax=0;
513 n1=0;
514 n2=ncons-1;
515
516 for(i=0; (i<ncons); i++) {
517 j=3*i;
518 a1=iatom[j+1];
519 a2=iatom[j+2];
520 nr=0;
521 for(k=0; (k<ncons); k++) {
522 b1=iatom[3*k+1];
523 b2=iatom[3*k+2];
524 if ((a1==b1 || a1==b2) || (a2==b1 || a2==b2))
525 if (i != k) nr++;
526 }
527 if (nr > *cmax) *cmax=nr;
528 type=iatom[j];
529 len =idef->iparams[type].shake.dA;
530 len1=idef->iparams[type].shake.dB;
531 if (nr <=4) {
532 (*bla1)[n1]=a1;
533 (*bla2)[n1]=a2;
534 (*bllen)[n1]=len;
535 (*bllen0)[n1]=len;
536 (*ddist)[n1]=len1-len;
537 n1++;
538 }
539 else {
540 (*bla1)[n2]=a1;
541 (*bla2)[n2]=a2;
542 (*bllen)[n2]=len;
543 (*bllen0)[n2]=len;
544 (*ddist)[n2]=len1-len;
545 n2--;
546 }
547 }
548
549 *ncm=ncons-n1;
550 nc4=(*cmax-4)*n1;
551
552 i=4*n1+(*cmax)*(*ncm);
553 snew(*blbnb,i);
554 snew(*blcc,i);
555 snew(*blm,i);
556
557 for(i=0; (i<ncons); i++) {
558 a1=(*bla1)[i];
559 a2=(*bla2)[i];
560 im1=md->invmass[a1];
561 im2=md->invmass[a2];
562 (*blc)[i]=invsqrt(im1+im2);
563 /* printf("%d %d %f\n",a1,a2,blist[i].len); */
564 nr=0;
565 /* printf("\n%d",i); */
566 for(k=0; (k<ncons); k++) {
567 b1=(*bla1)[k];
568 b2=(*bla2)[k];
569 if ((a1==b1 || a1==b2) || (a2==b1 || a2==b2))
570 if (i != k) {
571 if (i<n1)
572 (*blbnb)[4*i+nr]=k;
573 else
574 (*blbnb)[(*cmax)*i-nc4+nr]=k;
575 nr++;
576 /* printf(" %d",k); */
577 }
578 }
579 (*blnr)[i]=nr;
580 /* printf("\n"); */
581 }
582 fprintf(stdlog,"\nInitializing LINear Constraint Solver\n");
583 fprintf(stdlog,"%d constraints\nof which %d with more than 4 neighbours\n",ncons,*ncm);
584 fprintf(stdlog,"maximum number of bonds coupled to one bond is %d\n\n",*cmax);
585 fflush(stdlog);
586
587 for(b=0; (b<ncons); b++) {
588 i=(*bla1)[b];
589 j=(*bla2)[b];
590 nr=(*blnr)[b];
591 if (nr)
592 for(n=0; (n<nr);n++) {
593 if (b < n1)
594 k=(*blbnb)[4*b+n];
595 else
596 k=(*blbnb)[(*cmax)*b-nc4+n];
597 if (i==(*bla1)[k] || j==(*bla2)[k])
598 sign=-1;
599 else
600 sign=1;
601 if (i==(*bla1)[k] || i==(*bla2)[k])
602 cen=i;
603 else
604 cen=j;
605 if (ir->eI==eiMD) {
606 imcen=md->invmass[cen];
607 len=sign*imcen*(*blc)[b]*(*blc)[k];
608 }
609 if (ir->eI==eiLD)
610 len=sign*0.5;
611 if (b<n1)
612 (*blcc)[4*b+n]=len;
613 else
614 (*blcc)[(*cmax)*b-nc4+n]=len;
615 }
616 }
617
618 }
619 }
620
621 void dump_confs(int step,t_atoms *atoms,rvec x[],rvec xprime[],matrix box)
622 {
623 char buf[256];
624
625 sprintf(buf,"step%d.pdb",step-1);
626 write_sto_conf(buf,"one step before crash",atoms,x,NULL,box);
627 sprintf(buf,"step%d.pdb",step);
628 write_sto_conf(buf,"crashed",atoms,xprime,NULL,box);
629 fprintf(stdlog,"Wrote pdb files with previous and current coordinates\n");
630 fprintf(stderr,"Wrote pdb files with previous and current coordinates\n");
631 }
632
633 void update(int natoms, /* number of atoms in simulation */
634 int start,
635 int homenr, /* number of home particles */
636 int step,
637 real lambda,
638 real *dvdlambda, /* FEP stuff */
639 t_inputrec *ir, /* input record with constants */
640 bool bFirstStep,
641 t_mdatoms *md,
642 rvec x[], /* coordinates of home particles */
643 t_graph *graph,
644 rvec shift_vec[],
645 rvec force[], /* forces on home particles */
646 rvec delta_f[],
647 rvec vold[], /* Old velocities */
648 rvec v[], /* velocities of home particles */
649 rvec vt[], /* velocity at time t */
650 tensor pressure, /* instantaneous pressure tensor */
651 tensor box, /* instantaneous box lengths */
652 t_topology *top,
653 t_groups *grps,
654 tensor vir_part,
655 t_commrec *cr,
656 t_nrnb *nrnb,
657 bool bTYZ,
658 bool bDoUpdate,
659 t_edsamyn *edyn)
660 {
661 static char buf[256];
662 static bool bFirst=TRUE;
663 static int nblocks=0;
664 static int *sblock=NULL;
665 static int nsettle,settle_type;
666 static int *owptr;
667 static rvec *xprime,*x_unc=NULL;
668 static int ngtc,ngacc,ngfrz;
669 static rvec *lamb,*freezefac;
670 static int *bla1,*bla2,*blnr,*blbnb;
671 static rvec *r;
672 static real *bllen,*blc,*blcc,*blm,*tmp1,*tmp2,*tmp3,*lincslam,
673 *bllen0,*ddist;
674 static int ncm,cmax;
675 static t_edpar edpar;
676
677 t_idef *idef=&(top->idef);
678 double dt;
679 real dt_1,dt_2;
680 int i,n,m,g;
681 int ncons=0,nc;
682 int warn,p_imax,error;
683 real wang,p_max,p_rms;
684
685 set_led(UPDATE_LED);
686
687 if (bFirst) {
688 init_update(stdlog,top,ir,md,
689 start,homenr,
690 &nblocks,&sblock,
691 &nsettle,&owptr,&settle_type);
692
693 if (idef->il[F_SHAKE].nr) {
694 if (ir->eConstrAlg == estLINCS) {
695 please_cite(stdlog,"Hess97a");
696 init_lincs(stdlog,top,ir,md,
697 start,homenr,
698 &nblocks,&sblock,
699 &nsettle,&owptr,&settle_type,
700 &ncm,&cmax,
701 &r,&bla1,&bla2,&blnr,&blbnb,
702 &bllen,&blc,&blcc,&blm,&tmp1,&tmp2,&tmp3,&lincslam,
703 &bllen0,&ddist);
704 }
705 else
706 please_cite(stdlog,"Ryckaert77a");
707 }
708 if (edyn->bEdsam)
709 init_edsam(stdlog,top,md,start,homenr,&nblocks,&sblock,x,box,
710 edyn,&edpar);
711
712 /* Allocate memory for xold, original atomic positions
713 * and for xprime.
714 */
715 snew(xprime,natoms);
716 snew(x_unc,homenr);
717
718 /* Freeze Factor: If a dimension of a group has to be frozen,
719 * the corresponding freeze fac will be 0.0 otherwise 1.0
720 * This is implemented by multiplying the CHANGE in position
721 * by freeze fac (also in do_pcoupl)
722 *
723 * Coordinates in shake can be frozen by setting the invmass
724 * of a particle to 0.0 (===> Infinite mass!)
725 */
726 ngfrz=ir->opts.ngfrz;
727 snew(freezefac,ngfrz);
728 for(n=0; (n<ngfrz); n++)
729 for(m=0; (m<DIM); m++) {
730 freezefac[n][m]=(ir->opts.nFreeze[n][m]==0) ? 1.0 : 0.0;
731 /* printf("n %d m %d ff %g\n",n,m,freezefac[n][m]); */
732 }
733 /* for(i=0; (i<natoms); i++) */
734 /* printf("%d fg %d\n",i,md->cFREEZE[i]); */
735 /* Copy the pointer to the external acceleration in the opts */
736 ngacc=ir->opts.ngacc;
737
738 ngtc=ir->opts.ngtc;
739 snew(lamb,ir->opts.ngtc);
740
741 /* done with initializing */
742 bFirst=FALSE;
743 }
744
745 dt = ir->delta_t;
746 dt_1 = 1.0/dt;
747 dt_2 = 1.0/(dt*dt);
748
749 for(i=0; (i<ngtc); i++) {
750 real l=grps->tcstat[i].lambda;
751
752 if (bTYZ)
753 lamb[i][XX]=1;
754 else
755 lamb[i][XX]=l;
756 lamb[i][YY]=l;
757 lamb[i][ZZ]=l;
758 }
759
760 if (bDoUpdate) {
761 /* update mean velocities */
762 for (g=0; (g<ngacc); g++) {
763 copy_rvec(grps->grpstat[g].u,grps->grpstat[g].uold);
764 clear_rvec(grps->grpstat[g].u);
765 }
766
767 /* Now do the actual update of velocities and positions */
768 where();
769 dump_it_all(stdlog,"Before update",natoms,x,xprime,v,vold,force);
770 if (ir->eI==eiMD)
771 /* use normal version of update */
772 do_update(start,homenr,dt,
773 lamb,grps->grpstat,
774 ir->opts.acc,freezefac,
775 md->invmass,md->ptype,
776 md->cFREEZE,md->cACC,md->cTC,
777 x,xprime,v,vold,force);
778 else if (ir->eI==eiLD)
779 do_update_lang(start,homenr,dt,
780 freezefac,md->ptype,md->cFREEZE,
781 x,xprime,v,vold,force,
782 ir->ld_temp,ir->ld_fric,&ir->ld_seed);
783 else
784 fatal_error(0,"Don't know how to update coordinates");
785
786 where();
787 inc_nrnb(nrnb,eNR_UPDATE,homenr);
788 dump_it_all(stdlog,"After update",natoms,x,xprime,v,vold,force);
789 }
790 else {
791 /* If we're not updating we're doing shakefirst!
792 * In this case the extra coordinates are passed in v array
793 */
794 for(n=start; (n<start+homenr); n++) {
795 copy_rvec(v[n],xprime[n]);
796 }
797 }
798
799 /*
800 * Steps (7C, 8C)
801 * APPLY CONSTRAINTS:
802 * BLOCK SHAKE
803 */
804
805 if ((nblocks > 0) || (nsettle > 0)) {
806 /* Copy Unconstrained X to temp array */
807 for(n=start; (n<start+homenr); n++)
808 copy_rvec(xprime[n],x_unc[n-start]);
809
810
811 if (nblocks > 0) {
812 where();
813
814 nc=idef->il[F_SHAKE].nr/3;
815
816 if (ir->eConstrAlg == estSHAKE)
817 ncons=bshakef(stdlog,natoms,md->invmass,nblocks,sblock,idef,
818 ir,box,x,xprime,nrnb);
819
820 if (ir->eConstrAlg == estLINCS) {
821
822 if (ir->bPert) {
823 for(i=0;i<nc;i++)
824 bllen[i]=bllen0[i]+lambda*ddist[i];
825 }
826
827 wang=ir->LincsWarnAngle;
828
829 if (do_per_step(step,ir->nstLincsout))
830 cconerr(&p_max,&p_rms,&p_imax,xprime,nc,bla1,bla2,bllen);
831
832 if (ir->eI==eiMD) {
833 #ifdef USEF77
834 flincs(x[0],xprime[0],&nc,&ncm,&cmax,bla1,bla2,blnr,blbnb,
835 bllen,blc,blcc,blm,&ir->nProjOrder,
836 md->invmass,r[0],tmp1,tmp2,tmp3,&wang,&warn,
837 lincslam);
838 #else
839 clincs(x,xprime,nc,ncm,cmax,bla1,bla2,blnr,blbnb,
840 bllen,blc,blcc,blm,ir->nProjOrder,
841 md->invmass,r,tmp1,tmp2,tmp3,wang,&warn,lincslam);
842 #endif
843 if (ir->bPert) {
844 real dvdl=0;
845
846 for(i=0; (i<nc); i++)
847 dvdl+=lincslam[i]*dt_2*ddist[i];
848 *dvdlambda+=dvdl;
849 }
850 }
851
852 if (ir->eI==eiLD) {
853 #ifdef USEF77
854 flincsld(x[0],xprime[0],&nc,&ncm,&cmax,bla1,bla2,blnr,
855 blbnb,bllen,blcc,blm,&ir->nProjOrder,
856 r[0],tmp1,tmp2,tmp3,&wang,&warn);
857 #else
858 clincsld(x,xprime,nc,ncm,cmax,bla1,bla2,blnr,
859 blbnb,bllen,blcc,blm,ir->nProjOrder,
860 r,tmp1,tmp2,tmp3,wang,&warn);
861 #endif
862 }
863
864 if (do_per_step(step,ir->nstLincsout)) {
865 fprintf(stdlog,"Step %d\nRel. Constraint Deviation: Max between atoms RMS\n",step);
866 fprintf(stdlog," Before LINCS %.6f %6d %6d %.6f\n",
867 p_max,bla1[p_imax]+1,bla2[p_imax]+1,p_rms);
868 cconerr(&p_max,&p_rms,&p_imax,xprime,nc,bla1,bla2,bllen);
869 fprintf(stdlog," After LINCS %.6f %6d %6d %.6f\n\n",
870 p_max,bla1[p_imax]+1,bla2[p_imax]+1,p_rms);
871 }
872
873 if (warn > 0) {
874 cconerr(&p_max,&p_rms,&p_imax,xprime,nc,bla1,bla2,bllen);
875 sprintf(buf,"\nStep %d, time %g (ps) LINCS WARNING\n"
876 "relative constraint deviation after LINCS:\n"
877 "max %.6f (between atoms %d and %d) rms %.6f\n",
878 step,ir->init_t+step*dt,
879 p_max,bla1[p_imax]+1,bla2[p_imax]+1,p_rms);
880 fprintf(stdlog,"%s",buf);
881 fprintf(stderr,"%s",buf);
882 lincs_warning(x,xprime,nc,bla1,bla2,bllen,wang);
883 if (p_max > 0.5) {
884 dump_confs(step,&(top->atoms),x,xprime,box);
885 fatal_error(0,"Bond deviates more than half its own length");
886 }
887 }
888 }
889
890 where();
891
892
893 if (ncons == -1) {
894 dump_confs(step,&(top->atoms),x,xprime,box);
895 fprintf(stdlog,"SHAKE ERROR at step %d\n",step);
896 fatal_error(0,"SHAKE ERROR at step %d\n",step);
897 }
898
899 dump_it_all(stdlog,"After Shake",natoms,x,xprime,v,vold,force);
900 }
901
902 /* apply Essential Dynamics constraints when required */
903 if (edyn->bEdsam)
904 do_edsam(stdlog,top,ir,step,md,start,homenr,&nblocks,&sblock,xprime,x,
905 x_unc,force,box,edyn,&edpar,bDoUpdate);
906
907 if (nsettle > 0) {
908 int ow1;
909 real mO,mH,dOH,dHH;
910
911 ow1 = owptr[0];
912 mO = md->massA[ow1];
913 mH = md->massA[ow1+1];
914 dOH = top->idef.iparams[settle_type].settle.doh;
915 dHH = top->idef.iparams[settle_type].settle.dhh;
916 #ifdef USEF77
917 fsettle(&nsettle,owptr,x[0],xprime[0],&dOH,&dHH,&mO,&mH,&error);
918 #else
919 csettle(stdlog,nsettle,owptr,x[0],xprime[0],dOH,dHH,mO,mH,&error);
920 #endif
921 inc_nrnb(nrnb,eNR_SETTLE,nsettle);
922 if (error>=0) {
923 dump_confs(step,&(top->atoms),x,xprime,box);
924 sprintf(buf,"\nMolecule starting at atomnr. %d can not be settled, "
925 "step %d, time %g (ps)",
926 owptr[error]+1,step,ir->init_t+step*dt);
927 fprintf(stdlog,buf);
928 fatal_error(0,buf);
929 }
930 where();
931 }
932 if (bDoUpdate) {
933 real mdt_2;
934
935 for(n=start; (n<start+homenr); n++) {
936 mdt_2 = dt_2*md->massT[n];
937 do_both(x[n],x_unc[n-start],xprime[n],delta_f[n],
938 v[n],mdt_2,dt_1);
939 }
940 where();
941
942 inc_nrnb(nrnb,eNR_SHAKE_V,homenr);
943 dump_it_all(stdlog,"After Shake-V",natoms,x,xprime,v,vold,force);
944 where();
945
946 /* Calculate virial due to shake (for this proc) */
947 calc_vir(stdlog,homenr,&(x[start]),&(delta_f[start]),vir_part,cr);
948 inc_nrnb(nrnb,eNR_SHAKE_VIR,homenr);
949 where();
950 }
951 }
952
953
954 /* We must always unshift here, also if we did not shake */
955 where();
956 if ((graph->nnodes > 0) && bDoUpdate) {
957 unshift_x(graph,shift_vec,x,xprime);
958 inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);
959 for(n=start; (n<graph->start); n++)
960 copy_rvec(xprime[n],x[n]);
961 for(n=graph->start+graph->nnodes; (n<start+homenr); n++)
962 copy_rvec(xprime[n],x[n]);
963 }
964 else {
965 for(n=start; (n<start+homenr); n++)
966 copy_rvec(xprime[n],x[n]);
967 }
968 dump_it_all(stdlog,"After unshift",natoms,x,xprime,v,vold,force);
969 where();
970
971 if (bDoUpdate) {
972 update_grps(start,homenr,grps,&(ir->opts),v,md);
973 do_pcoupl(ir,step,pressure,box,start,homenr,x,md->cFREEZE,nrnb,freezefac);
974 where();
975 }
976
977 clr_led(UPDATE_LED);
978 }
got merge a long time ago