?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed a bug in the pdb-writing code.
[gromacs.git] / src / kernel / do_gct.c
blob28618c7f5688704e063f5a7ad4c056199a5e8bf5
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 3.1
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
16 *
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
23 *
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
26 *
27 * For more info, check our website at http://www.gromacs.org
28 *
29 * And Hey:
30 * GROningen Mixture of Alchemy and Childrens' Stories
31 */
32 static char *SRCID_do_gct_c = "$Id$";
33 #include "typedefs.h"
34 #include "xmdrun.h"
35 #include "block_tx.h"
36 #include "futil.h"
37 #include "xvgr.h"
38 #include "macros.h"
39 #include "physics.h"
40 #include "network.h"
41 #include "smalloc.h"
42 #include "mdrun.h"
43 #include "string2.h"
44 #include "readinp.h"
45 #include "filenm.h"
46 #include "update.h"
47 #include "vec.h"
48 #include "main.h"
49 #include "txtdump.h"
50
51 /*#define DEBUGGCT*/
52
53 t_coupl_rec *init_coupling(FILE *log,int nfile,t_filenm fnm[],
54 t_commrec *cr,t_forcerec *fr,
55 t_mdatoms *md,t_idef *idef)
56 {
57 int i,nc,index,j;
58 int ati,atj;
59 t_coupl_rec *tcr;
60
61 snew(tcr,1);
62 if (MASTER(cr)) {
63 read_gct (opt2fn("-j",nfile,fnm), tcr);
64 write_gct(opt2fn("-jo",nfile,fnm),tcr,idef);
65 }
66 copy_ff(tcr,fr,md,idef);
67
68 /* Update all processors with coupling info */
69 if (PAR(cr))
70 comm_tcr(log,cr,&tcr);
71
72 return tcr;
73 }
74
75 static real Ecouple(t_coupl_rec *tcr,real ener[])
76 {
77 if (tcr->bInter)
78 return ener[F_SR]+ener[F_LJ]+ener[F_LR]+ener[F_LJLR];
79 else
80 return ener[F_EPOT];
81 }
82
83 static char *mk_gct_nm(char *fn,int ftp,int ati,int atj)
84 {
85 static char buf[256];
86
87 strcpy(buf,fn);
88 if (atj == -1)
89 sprintf(buf+strlen(fn)-4,"%d.%s",ati,ftp2ext(ftp));
90 else
91 sprintf(buf+strlen(fn)-4,"%d_%d.%s",ati,atj,ftp2ext(ftp));
92
93 return buf;
94 }
95
96 static void pr_ff(t_coupl_rec *tcr,real time,t_idef *idef,
97 t_commrec *cr,int nfile,t_filenm fnm[])
98 {
99 static FILE *prop;
100 static FILE **out=NULL;
101 static FILE **qq=NULL;
102 static FILE **ip=NULL;
103 t_coupl_LJ *tclj;
104 t_coupl_BU *tcbu;
105 char buf[256];
106 char *leg[] = { "C12", "C6" };
107 char *bleg[] = { "A", "B", "C" };
108 char **raleg;
109 int i,j,index;
110
111 if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL)) {
112 prop=xvgropen(opt2fn("-runav",nfile,fnm),
113 "Properties and Running Averages","Time (ps)","");
114 snew(raleg,2*eoObsNR);
115 for(i=j=0; (i<eoObsNR); i++) {
116 if (tcr->bObsUsed[i]) {
117 raleg[j++] = strdup(eoNames[i]);
118 sprintf(buf,"RA-%s",eoNames[i]);
119 raleg[j++] = strdup(buf);
120 }
121 }
122 xvgr_legend(prop,j,raleg);
123 for(i=0; (i<j); i++)
124 sfree(raleg[i]);
125 sfree(raleg);
126
127 if (tcr->nLJ) {
128 snew(out,tcr->nLJ);
129 for(i=0; (i<tcr->nLJ); i++) {
130 if (tcr->tcLJ[i].bPrint) {
131 tclj = &(tcr->tcLJ[i]);
132 out[i] =
133 xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),
134 efXVG,tclj->at_i,tclj->at_j),
135 "General Coupling Lennard Jones","Time (ps)",
136 "Force constant (units)");
137 fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
138 tclj->at_i,tclj->at_j);
139 xvgr_legend(out[i],asize(leg),leg);
140 fflush(out[i]);
141 }
142 }
143 }
144 else if (tcr->nBU) {
145 snew(out,tcr->nBU);
146 for(i=0; (i<tcr->nBU); i++) {
147 if (tcr->tcBU[i].bPrint) {
148 tcbu=&(tcr->tcBU[i]);
149 out[i] =
150 xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
151 tcbu->at_i,tcbu->at_j),
152 "General Coupling Buckingham","Time (ps)",
153 "Force constant (units)");
154 fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
155 tcbu->at_i,tcbu->at_j);
156 xvgr_legend(out[i],asize(bleg),bleg);
157 fflush(out[i]);
158 }
159 }
160 }
161 snew(qq,tcr->nQ);
162 for(i=0; (i<tcr->nQ); i++) {
163 if (tcr->tcQ[i].bPrint) {
164 qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
165 tcr->tcQ[i].at_i,-1),
166 "General Coupling Charge","Time (ps)","Charge (e)");
167 fprintf(qq[i],"@ subtitle \"Type %d\"\n",tcr->tcQ[i].at_i);
168 fflush(qq[i]);
169 }
170 }
171 snew(ip,tcr->nIP);
172 for(i=0; (i<tcr->nIP); i++) {
173 sprintf(buf,"gctIP%d",tcr->tIP[i].type);
174 ip[i]=xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,0,-1),
175 "General Coupling iparams","Time (ps)","ip ()");
176 index=tcr->tIP[i].type;
177 fprintf(ip[i],"@ subtitle \"Coupling to %s\"\n",
178 interaction_function[idef->functype[index]].longname);
179 fflush(ip[i]);
180 }
181 }
182 /* Write properties to file */
183 fprintf(prop,"%10.3f",time);
184 for(i=0; (i<eoObsNR); i++)
185 if (tcr->bObsUsed[i])
186 fprintf(prop," %10.3e %10.3e",tcr->act_value[i],tcr->av_value[i]);
187 fprintf(prop,"\n");
188 fflush(prop);
189
190 for(i=0; (i<tcr->nLJ); i++) {
191 tclj=&(tcr->tcLJ[i]);
192 if (tclj->bPrint) {
193 fprintf(out[i],"%14.7e %14.7e %14.7e\n",
194 time,tclj->c12,tclj->c6);
195 fflush(out[i]);
196 }
197 }
198 for(i=0; (i<tcr->nBU); i++) {
199 tcbu=&(tcr->tcBU[i]);
200 if (tcbu->bPrint) {
201 fprintf(out[i],"%14.7e %14.7e %14.7e %14.7e\n",
202 time,tcbu->a,tcbu->b,tcbu->c);
203 fflush(out[i]);
204 }
205 }
206 for(i=0; (i<tcr->nQ); i++) {
207 if (tcr->tcQ[i].bPrint) {
208 fprintf(qq[i],"%14.7e %14.7e\n",time,tcr->tcQ[i].Q);
209 fflush(qq[i]);
210 }
211 }
212 for(i=0; (i<tcr->nIP); i++) {
213 fprintf(ip[i],"%10g ",time);
214 index=tcr->tIP[i].type;
215 switch(idef->functype[index]) {
216 case F_BONDS:
217 fprintf(ip[i],"%10g %10g\n",tcr->tIP[i].iprint.harmonic.krA,
218 tcr->tIP[i].iprint.harmonic.rA);
219 break;
220 default:
221 break;
222 }
223 fflush(ip[i]);
224 }
225 }
226
227 static void pr_dev(t_coupl_rec *tcr,
228 real t,real dev[eoObsNR],t_commrec *cr,int nfile,t_filenm fnm[])
229 {
230 static FILE *fp=NULL;
231 char **ptr;
232 int i,j;
233
234 if (!fp) {
235 fp=xvgropen(opt2fn("-devout",nfile,fnm),
236 "Deviations from target value","Time (ps)","");
237 snew(ptr,eoObsNR);
238 for(i=j=0; (i<eoObsNR); i++)
239 if (tcr->bObsUsed[i])
240 ptr[j++] = eoNames[i];
241 xvgr_legend(fp,j,ptr);
242 sfree(ptr);
243 }
244 fprintf(fp,"%10.3f",t);
245 for(i=0; (i<eoObsNR); i++)
246 if (tcr->bObsUsed[i])
247 fprintf(fp," %10.3e",dev[i]);
248 fprintf(fp,"\n");
249 fflush(fp);
250 }
251
252 static void upd_nbfplj(FILE *log,real *nbfp,int atnr,real f6[],real f12[])
253 {
254 int n,m,k;
255
256 /* Update the nonbonded force parameters */
257 for(k=n=0; (n<atnr); n++) {
258 for(m=0; (m<atnr); m++,k++) {
259 C6 (nbfp,atnr,n,m) *= f6[k];
260 C12(nbfp,atnr,n,m) *= f12[k];
261 }
262 }
263 }
264
265 static void upd_nbfpbu(FILE *log,real *nbfp,int atnr,
266 real fa[],real fb[],real fc[])
267 {
268 int n,m,k;
269
270 /* Update the nonbonded force parameters */
271 for(k=n=0; (n<atnr); n++) {
272 for(m=0; (m<atnr); m++,k++) {
273 BHAMA(nbfp,atnr,n,m) *= fa[k];
274 BHAMB(nbfp,atnr,n,m) *= fb[k];
275 BHAMC(nbfp,atnr,n,m) *= fc[k];
276 }
277 }
278 }
279
280 void gprod(t_commrec *cr,int n,real f[])
281 {
282 /* Compute the global product of all elements in an array
283 * such that after gprod f[i] = PROD_j=1,nprocs f[i][j]
284 */
285 static int nbuf=0;
286 static real *buf[2] = { NULL, NULL };
287 int i,j,cur=0;
288 #define next (1-cur)
289
290 if (n > nbuf) {
291 nbuf = n;
292 srenew(buf[cur], nbuf);
293 srenew(buf[next],nbuf);
294 }
295
296 for(j=0; (j<n); j++)
297 buf[cur][j] = f[j];
298
299 for(i=0; (i<cr->nnodes-1); i++) {
300 gmx_tx(cr->left, array(buf[cur],n));
301 gmx_rx(cr->right,array(buf[next],n));
302 gmx_wait(cr->left,cr->right);
303 /* Multiply f by factor read */
304 for(j=0; (j<n); j++)
305 f[j] *= buf[next][j];
306 /* Swap buffers */
307 cur = next;
308 }
309 #undef next
310 }
311
312 static void set_factor_matrix(int ntypes,real f[],real fmult,int ati,int atj)
313 {
314 #define FMIN 0.95
315 #define FMAX 1.05
316 int i;
317
318 fmult = min(FMAX,max(FMIN,fmult));
319 if (atj != -1) {
320 f[ntypes*ati+atj] *= fmult;
321 f[ntypes*atj+ati] *= fmult;
322 }
323 else {
324 for(i=0; (i<ntypes); i++) {
325 f[ntypes*ati+i] *= fmult;
326 f[ntypes*i+ati] *= fmult;
327 }
328 }
329 #undef FMIN
330 #undef FMAX
331 }
332
333 static real calc_deviation(real xav,real xt,real x0)
334 {
335 /* This may prevent overshooting in GCT coupling... */
336
337 /* real dev;
338
339 if (xav > x0) {
340 if (xt > x0)
341 dev = min(xav-x0,xt-x0);
342 else
343 dev = 0;
344 }
345 else {
346 if (xt < x0)
347 dev = max(xav-x0,xt-x0);
348 else
349 dev = 0;
350 }
351 */
352 return x0-xav;
353 }
354
355 static real calc_dist(FILE *log,rvec x[])
356 {
357 static bool bFirst=TRUE;
358 static bool bDist;
359 static int i1,i2;
360 char *buf;
361 rvec dx;
362
363 if (bFirst) {
364 if ((buf = getenv("DISTGCT")) == NULL)
365 bDist = FALSE;
366 else {
367 bDist = (sscanf(buf,"%d%d",&i1,&i2) == 2);
368 if (bDist)
369 fprintf(log,"GCT: Will couple to distance between %d and %d\n",i1,i2);
370 else
371 fprintf(log,"GCT: Will not couple to distances\n");
372 }
373 bFirst = FALSE;
374 }
375 if (bDist) {
376 rvec_sub(x[i1],x[i2],dx);
377 return norm(dx);
378 }
379 else
380 return 0.0;
381 }
382
383 real run_aver(real old,real cur,int step,int nmem)
384 {
385 nmem = max(1,nmem);
386
387 return ((nmem-1)*old+cur)/nmem;
388 }
389
390 static void set_act_value(t_coupl_rec *tcr,int index,real val,int step)
391 {
392 tcr->act_value[index] = val;
393 tcr->av_value[index] = run_aver(tcr->av_value[index],val,step,tcr->nmemory);
394 }
395
396 static void upd_f_value(FILE *log,int atnr,real xi,real dt,real factor,
397 real ff[],int ati,int atj)
398 {
399 real fff;
400
401 if (xi != 0) {
402 fff = 1 + (dt/xi) * factor;
403 if (fff > 0)
404 set_factor_matrix(atnr,ff,sqrt(fff),ati,atj);
405 }
406 }
407
408 static void dump_fm(FILE *fp,int n,real f[],char *s)
409 {
410 int i,j;
411
412 fprintf(fp,"Factor matrix (all numbers -1) %s\n",s);
413 for(i=0; (i<n); i++) {
414 for(j=0; (j<n); j++)
415 fprintf(fp," %10.3e",f[n*i+j]-1.0);
416 fprintf(fp,"\n");
417 }
418 }
419
420 void do_coupling(FILE *log,int nfile,t_filenm fnm[],
421 t_coupl_rec *tcr,real t,int step,real ener[],
422 t_forcerec *fr,t_inputrec *ir,bool bMaster,
423 t_mdatoms *md,t_idef *idef,real mu_aver,int nmols,
424 t_commrec *cr,matrix box,tensor virial,
425 tensor pres,rvec mu_tot,
426 rvec x[],rvec f[],bool bDoIt)
427 {
428 #define enm2Debye 48.0321
429 #define d2e(x) (x)/enm2Debye
430 #define enm2kjmol(x) (x)*0.0143952 /* = 2.0*4.0*M_PI*EPSILON0 */
431
432 static real *f6,*f12,*fa,*fb,*fc,*fq;
433 static bool bFirst = TRUE;
434
435 int i,j,ati,atj,atnr2,type,ftype;
436 real deviation[eoObsNR],prdev[eoObsNR],epot0,dist,rmsf;
437 real ff6,ff12,ffa,ffb,ffc,ffq,factor,dt,mu_ind;
438 real Epol,Eintern,Virial,muabs,xiH=-1,xiS=-1,xi6,xi12;
439 rvec fmol[2];
440 bool bTest,bPrint;
441 t_coupl_LJ *tclj;
442 t_coupl_BU *tcbu;
443 t_coupl_Q *tcq;
444 t_coupl_iparams *tip;
445
446 atnr2 = idef->atnr * idef->atnr;
447 if (bFirst) {
448 if (PAR(cr))
449 fprintf(log,"GCT: this is parallel\n");
450 else
451 fprintf(log,"GCT: this is not parallel\n");
452 fflush(log);
453 snew(f6, atnr2);
454 snew(f12,atnr2);
455 snew(fa, atnr2);
456 snew(fb, atnr2);
457 snew(fc, atnr2);
458 snew(fq, idef->atnr);
459
460 if (tcr->bVirial) {
461 int nrdf = 0;
462 real TTT = 0;
463 real Vol = det(box);
464
465 for(i=0; (i<ir->opts.ngtc); i++) {
466 nrdf += ir->opts.nrdf[i];
467 TTT += ir->opts.nrdf[i]*ir->opts.ref_t[i];
468 }
469 TTT /= nrdf;
470
471 /* Calculate reference virial from reference temperature and pressure */
472 tcr->ref_value[eoVir] = 0.5*BOLTZ*nrdf*TTT - (3.0/2.0)*
473 Vol*tcr->ref_value[eoPres];
474
475 fprintf(log,"GCT: TTT = %g, nrdf = %d, vir0 = %g, Vol = %g\n",
476 TTT,nrdf,tcr->ref_value[eoVir],Vol);
477 fflush(log);
478 }
479 bFirst = FALSE;
480 }
481
482 bPrint = MASTER(cr) && do_per_step(step,ir->nstlog);
483 dt = ir->delta_t;
484
485 /* Initiate coupling to the reference pressure and temperature to start
486 * coupling slowly.
487 */
488 if (step == 0) {
489 for(i=0; (i<eoObsNR); i++)
490 tcr->av_value[i] = tcr->ref_value[i];
491 if ((tcr->ref_value[eoDipole]) != 0.0) {
492 mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
493 Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
494 tcr->av_value[eoEpot] -= Epol;
495 fprintf(log,"GCT: mu_aver = %g(D), mu_ind = %g(D), Epol = %g (kJ/mol)\n",
496 mu_aver*enm2Debye,mu_ind*enm2Debye,Epol);
497 }
498 }
499
500 /* We want to optimize the LJ params, usually to the Vaporization energy
501 * therefore we only count intermolecular degrees of freedom.
502 * Note that this is now optional. switch UseEinter to yes in your gct file
503 * if you want this.
504 */
505 dist = calc_dist(log,x);
506 muabs = norm(mu_tot);
507 Eintern = Ecouple(tcr,ener)/nmols;
508 Virial = virial[XX][XX]+virial[YY][YY]+virial[ZZ][ZZ];
509
510 /*calc_force(md->nr,f,fmol);*/
511 clear_rvec(fmol[0]);
512
513 /* Use a memory of tcr->nmemory steps, so we actually couple to the
514 * average observable over the last tcr->nmemory steps. This may help
515 * in avoiding local minima in parameter space.
516 */
517 set_act_value(tcr,eoPres, ener[F_PRES],step);
518 set_act_value(tcr,eoEpot, Eintern, step);
519 set_act_value(tcr,eoVir, Virial, step);
520 set_act_value(tcr,eoDist, dist, step);
521 set_act_value(tcr,eoMu, muabs, step);
522 set_act_value(tcr,eoFx, fmol[0][XX], step);
523 set_act_value(tcr,eoFy, fmol[0][YY], step);
524 set_act_value(tcr,eoFz, fmol[0][ZZ], step);
525 set_act_value(tcr,eoPx, pres[XX][XX],step);
526 set_act_value(tcr,eoPy, pres[YY][YY],step);
527 set_act_value(tcr,eoPz, pres[ZZ][ZZ],step);
528
529 epot0 = tcr->ref_value[eoEpot];
530 /* If dipole != 0.0 assume we want to use polarization corrected coupling */
531 if ((tcr->ref_value[eoDipole]) != 0.0) {
532 mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
533
534 Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
535
536 epot0 -= Epol;
537 if (debug) {
538 fprintf(debug,"mu_ind: %g (%g D) mu_aver: %g (%g D)\n",
539 mu_ind,mu_ind*enm2Debye,mu_aver,mu_aver*enm2Debye);
540 fprintf(debug,"Eref %g Epol %g Erunav %g Eact %g\n",
541 tcr->ref_value[eoEpot],Epol,tcr->av_value[eoEpot],
542 tcr->act_value[eoEpot]);
543 }
544 }
545
546 if (bPrint) {
547 pr_ff(tcr,t,idef,cr,nfile,fnm);
548 }
549 /* Calculate the deviation of average value from the target value */
550 for(i=0; (i<eoObsNR); i++) {
551 deviation[i] = calc_deviation(tcr->av_value[i],tcr->act_value[i],
552 tcr->ref_value[i]);
553 prdev[i] = tcr->ref_value[i] - tcr->act_value[i];
554 }
555 deviation[eoEpot] = calc_deviation(tcr->av_value[eoEpot],tcr->act_value[eoEpot],
556 epot0);
557 prdev[eoEpot] = epot0 - tcr->act_value[eoEpot];
558
559 if (bPrint)
560 pr_dev(tcr,t,prdev,cr,nfile,fnm);
561
562 /* First set all factors to 1 */
563 for(i=0; (i<atnr2); i++) {
564 f6[i] = f12[i] = fa[i] = fb[i] = fc[i] = 1.0;
565 }
566 for(i=0; (i<idef->atnr); i++)
567 fq[i] = 1.0;
568
569 /* Now compute the actual coupling compononents */
570 if (!fr->bBHAM) {
571 if (bDoIt) {
572 for(i=0; (i<tcr->nLJ); i++) {
573 tclj=&(tcr->tcLJ[i]);
574
575 ati=tclj->at_i;
576 atj=tclj->at_j;
577
578 ff6 = ff12 = 1.0;
579
580 if (tclj->eObs == eoForce) {
581 fatal_error(0,"Hack code for this to work again ");
582 if (debug)
583 fprintf(debug,"Have computed derivatives: xiH = %g, xiS = %g\n",xiH,xiS);
584 if (ati == 1) {
585 /* Hydrogen */
586 ff12 += xiH;
587 }
588 else if (ati == 2) {
589 /* Shell */
590 ff12 += xiS;
591 }
592 else
593 fatal_error(0,"No H, no Shell, edit code at %s, line %d\n",
594 __FILE__,__LINE__);
595 if (ff6 > 0)
596 set_factor_matrix(idef->atnr,f6, sqrt(ff6), ati,atj);
597 if (ff12 > 0)
598 set_factor_matrix(idef->atnr,f12,sqrt(ff12),ati,atj);
599 }
600 else {
601 if (debug)
602 fprintf(debug,"tcr->tcLJ[%d].xi_6 = %g, xi_12 = %g deviation = %g\n",i,
603 tclj->xi_6,tclj->xi_12,deviation[tclj->eObs]);
604 factor=deviation[tclj->eObs];
605
606 upd_f_value(log,idef->atnr,tclj->xi_6, dt,factor,f6, ati,atj);
607 upd_f_value(log,idef->atnr,tclj->xi_12,dt,factor,f12,ati,atj);
608 }
609 }
610 }
611 if (PAR(cr)) {
612 gprod(cr,atnr2,f6);
613 gprod(cr,atnr2,f12);
614 #ifdef DEBUGGCT
615 dump_fm(log,idef->atnr,f6,"f6");
616 dump_fm(log,idef->atnr,f12,"f12");
617 #endif
618 }
619 upd_nbfplj(log,fr->nbfp,idef->atnr,f6,f12);
620
621 /* Copy for printing */
622 for(i=0; (i<tcr->nLJ); i++) {
623 tclj=&(tcr->tcLJ[i]);
624 ati = tclj->at_i;
625 atj = tclj->at_j;
626 if (atj == -1)
627 atj = ati;
628 tclj->c6 = C6(fr->nbfp,fr->ntype,ati,atj);
629 tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj);
630 }
631 }
632 else {
633 if (bDoIt) {
634 for(i=0; (i<tcr->nBU); i++) {
635 tcbu = &(tcr->tcBU[i]);
636 factor = deviation[tcbu->eObs];
637 ati = tcbu->at_i;
638 atj = tcbu->at_j;
639
640 upd_f_value(log,idef->atnr,tcbu->xi_a,dt,factor,fa,ati,atj);
641 upd_f_value(log,idef->atnr,tcbu->xi_b,dt,factor,fb,ati,atj);
642 upd_f_value(log,idef->atnr,tcbu->xi_c,dt,factor,fc,ati,atj);
643 }
644 }
645 if (PAR(cr)) {
646 gprod(cr,atnr2,fa);
647 gprod(cr,atnr2,fb);
648 gprod(cr,atnr2,fc);
649 }
650 upd_nbfpbu(log,fr->nbfp,idef->atnr,fa,fb,fc);
651 /* Copy for printing */
652 for(i=0; (i<tcr->nBU); i++) {
653 tcbu=&(tcr->tcBU[i]);
654 ati = tcbu->at_i;
655 atj = tcbu->at_j;
656 if (atj == -1)
657 atj = ati;
658 tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
659 tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
660 tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj);
661 if (debug)
662 fprintf(debug,"buck (type=%d) = %e, %e, %e\n",
663 tcbu->at_i,tcbu->a,tcbu->b,tcbu->c);
664 }
665 }
666 if (bDoIt) {
667 for(i=0; (i<tcr->nQ); i++) {
668 tcq=&(tcr->tcQ[i]);
669 if (tcq->xi_Q)
670 ffq = 1.0 + (dt/tcq->xi_Q) * deviation[tcq->eObs];
671 else
672 ffq = 1.0;
673 fq[tcq->at_i] *= ffq;
674 }
675 }
676 if (PAR(cr))
677 gprod(cr,idef->atnr,fq);
678
679 for(j=0; (j<md->nr); j++) {
680 md->chargeA[j] *= fq[md->typeA[j]];
681 }
682 for(i=0; (i<tcr->nQ); i++) {
683 tcq=&(tcr->tcQ[i]);
684 for(j=0; (j<md->nr); j++) {
685 if (md->typeA[j] == tcq->at_i) {
686 tcq->Q = md->chargeA[j];
687 break;
688 }
689 }
690 if (j == md->nr)
691 fatal_error(0,"Coupling type %d not found",tcq->at_i);
692 }
693 for(i=0; (i<tcr->nIP); i++) {
694 tip = &(tcr->tIP[i]);
695 type = tip->type;
696 ftype = idef->functype[type];
697 factor = dt*deviation[tip->eObs];
698
699 switch(ftype) {
700 case F_BONDS:
701 if (tip->xi.harmonic.krA) idef->iparams[type].harmonic.krA *= (1+factor/tip->xi.harmonic.krA);
702 if (tip->xi.harmonic.rA) idef->iparams[type].harmonic.rA *= (1+factor/tip->xi.harmonic.rA);
703 break;
704 default:
705 break;
706 }
707 tip->iprint=idef->iparams[type];
708 }
709 }
710
The ultimate molecular dynamics simulation package