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