2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
60 t_coupl_rec *init_coupling(FILE *log,int nfile, const t_filenm fnm[],
61 t_commrec *cr,t_forcerec *fr,
62 t_mdatoms *md,t_idef *idef)
70 read_gct (opt2fn("-j",nfile,fnm), tcr);
71 write_gct(opt2fn("-jo",nfile,fnm),tcr,idef);
73 /* Update all processors with coupling info */
75 comm_tcr(log,cr,&tcr);
77 /* Copy information from the coupling to the force field stuff */
78 copy_ff(tcr,fr,md,idef);
83 static real Ecouple(t_coupl_rec *tcr,real ener[])
86 return ener[F_COUL_SR]+ener[F_LJ]+ener[F_COUL_LR]+ener[F_LJ_LR]+
87 ener[F_RF_EXCL]+ener[F_COUL_RECIP];
92 static char *mk_gct_nm(const char *fn,int ftp,int ati,int atj)
98 sprintf(buf+strlen(fn)-4,"%d.%s",ati,ftp2ext(ftp));
100 sprintf(buf+strlen(fn)-4,"%d_%d.%s",ati,atj,ftp2ext(ftp));
105 static void pr_ff(t_coupl_rec *tcr,real time,t_idef *idef,
106 t_commrec *cr,int nfile,const t_filenm fnm[],
107 const output_env_t oenv)
110 static FILE **out=NULL;
111 static FILE **qq=NULL;
112 static FILE **ip=NULL;
116 const char *leg[] = { "C12", "C6" };
117 const char *eleg[] = { "Epsilon", "Sigma" };
118 const char *bleg[] = { "A", "B", "C" };
122 if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL)) {
123 prop=xvgropen(opt2fn("-runav",nfile,fnm),
124 "Properties and Running Averages","Time (ps)","",oenv);
125 snew(raleg,2*eoObsNR);
126 for(i=j=0; (i<eoObsNR); i++) {
127 if (tcr->bObsUsed[i]) {
128 raleg[j++] = strdup(eoNames[i]);
129 sprintf(buf,"RA-%s",eoNames[i]);
130 raleg[j++] = strdup(buf);
133 xvgr_legend(prop,j,(const char**)raleg,oenv);
140 for(i=0; (i<tcr->nLJ); i++) {
141 if (tcr->tcLJ[i].bPrint) {
142 tclj = &(tcr->tcLJ[i]);
144 xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),
145 efXVG,tclj->at_i,tclj->at_j),
146 "General Coupling Lennard Jones","Time (ps)",
147 "Force constant (units)",oenv);
148 fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
149 tclj->at_i,tclj->at_j);
150 if (tcr->combrule == 1)
151 xvgr_legend(out[i],asize(leg),leg,oenv);
153 xvgr_legend(out[i],asize(eleg),eleg,oenv);
160 for(i=0; (i<tcr->nBU); i++) {
161 if (tcr->tcBU[i].bPrint) {
162 tcbu=&(tcr->tcBU[i]);
164 xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
165 tcbu->at_i,tcbu->at_j),
166 "General Coupling Buckingham","Time (ps)",
167 "Force constant (units)",oenv);
168 fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
169 tcbu->at_i,tcbu->at_j);
170 xvgr_legend(out[i],asize(bleg),bleg,oenv);
176 for(i=0; (i<tcr->nQ); i++) {
177 if (tcr->tcQ[i].bPrint) {
178 qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
179 tcr->tcQ[i].at_i,-1),
180 "General Coupling Charge","Time (ps)","Charge (e)",
182 fprintf(qq[i],"@ subtitle \"Type %d\"\n",tcr->tcQ[i].at_i);
187 for(i=0; (i<tcr->nIP); i++) {
188 sprintf(buf,"gctIP%d",tcr->tIP[i].type);
189 ip[i]=xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,0,-1),
190 "General Coupling iparams","Time (ps)","ip ()",oenv);
191 index=tcr->tIP[i].type;
192 fprintf(ip[i],"@ subtitle \"Coupling to %s\"\n",
193 interaction_function[idef->functype[index]].longname);
197 /* Write properties to file */
198 fprintf(prop,"%10.3f",time);
199 for(i=0; (i<eoObsNR); i++)
200 if (tcr->bObsUsed[i])
201 fprintf(prop," %10.3e %10.3e",tcr->act_value[i],tcr->av_value[i]);
205 for(i=0; (i<tcr->nLJ); i++) {
206 tclj=&(tcr->tcLJ[i]);
208 if (tcr->combrule == 1)
209 fprintf(out[i],"%14.7e %14.7e %14.7e\n",
210 time,tclj->c12,tclj->c6);
212 double sigma = pow(tclj->c12/tclj->c6,1.0/6.0);
213 double epsilon = 0.25*sqr(tclj->c6)/tclj->c12;
214 fprintf(out[i],"%14.7e %14.7e %14.7e\n",
220 for(i=0; (i<tcr->nBU); i++) {
221 tcbu=&(tcr->tcBU[i]);
223 fprintf(out[i],"%14.7e %14.7e %14.7e %14.7e\n",
224 time,tcbu->a,tcbu->b,tcbu->c);
228 for(i=0; (i<tcr->nQ); i++) {
229 if (tcr->tcQ[i].bPrint) {
230 fprintf(qq[i],"%14.7e %14.7e\n",time,tcr->tcQ[i].Q);
234 for(i=0; (i<tcr->nIP); i++) {
235 fprintf(ip[i],"%10g ",time);
236 index=tcr->tIP[i].type;
237 switch(idef->functype[index]) {
239 fprintf(ip[i],"%10g %10g\n",tcr->tIP[i].iprint.harmonic.krA,
240 tcr->tIP[i].iprint.harmonic.rA);
249 static void pr_dev(t_coupl_rec *tcr,
250 real t,real dev[eoObsNR],t_commrec *cr,int nfile,
251 const t_filenm fnm[],const output_env_t oenv)
253 static FILE *fp=NULL;
258 fp=xvgropen(opt2fn("-devout",nfile,fnm),
259 "Deviations from target value","Time (ps)","",oenv);
261 for(i=j=0; (i<eoObsNR); i++)
262 if (tcr->bObsUsed[i])
263 ptr[j++] = strdup(eoNames[i]);
264 xvgr_legend(fp,j,(const char**)ptr,oenv);
269 fprintf(fp,"%10.3f",t);
270 for(i=0; (i<eoObsNR); i++)
271 if (tcr->bObsUsed[i])
272 fprintf(fp," %10.3e",dev[i]);
277 static void upd_nbfplj(FILE *log,real *nbfp,int atnr,real f6[],real f12[],
280 double *sigma,*epsilon,c6,c12,eps,sig,sig6;
283 /* Update the nonbonded force parameters */
286 for(k=n=0; (n<atnr); n++) {
287 for(m=0; (m<atnr); m++,k++) {
288 C6 (nbfp,atnr,n,m) *= f6[k];
289 C12(nbfp,atnr,n,m) *= f12[k];
295 /* Convert to sigma and epsilon */
298 for(n=0; (n<atnr); n++) {
300 c6 = C6 (nbfp,atnr,n,n) * f6[k];
301 c12 = C12(nbfp,atnr,n,n) * f12[k];
302 if ((c6 == 0) || (c12 == 0))
303 gmx_fatal(FARGS,"You can not use combination rule %d with zero C6 (%f) or C12 (%f)",combrule,c6,c12);
304 sigma[n] = pow(c12/c6,1.0/6.0);
305 epsilon[n] = 0.25*(c6*c6/c12);
307 for(k=n=0; (n<atnr); n++) {
308 for(m=0; (m<atnr); m++,k++) {
309 eps = sqrt(epsilon[n]*epsilon[m]);
311 sig = 0.5*(sigma[n]+sigma[m]);
313 sig = sqrt(sigma[n]*sigma[m]);
315 /* nbfp now includes the 6.0/12.0 derivative prefactors */
316 C6 (nbfp,atnr,n,m) = 4*eps*sig6/6.0;
317 C12(nbfp,atnr,n,m) = 4*eps*sig6*sig6/12.0;
324 gmx_fatal(FARGS,"Combination rule should be 1,2 or 3 instead of %d",
329 static void upd_nbfpbu(FILE *log,real *nbfp,int atnr,
330 real fa[],real fb[],real fc[])
334 /* Update the nonbonded force parameters */
335 for(k=n=0; (n<atnr); n++) {
336 for(m=0; (m<atnr); m++,k++) {
337 BHAMA(nbfp,atnr,n,m) *= fa[k];
338 BHAMB(nbfp,atnr,n,m) *= fb[k];
339 BHAMC(nbfp,atnr,n,m) *= fc[k];
344 void gprod(t_commrec *cr,int n,real f[])
346 /* Compute the global product of all elements in an array
347 * such that after gprod f[i] = PROD_j=1,nprocs f[i][j]
350 static real *buf=NULL;
360 MPI_Allreduce(f,buf,n,MPI_DOUBLE,MPI_PROD,cr->mpi_comm_mygroup);
362 MPI_Allreduce(f,buf,n,MPI_FLOAT, MPI_PROD,cr->mpi_comm_mygroup);
369 static void set_factor_matrix(int ntypes,real f[],real fmult,int ati,int atj)
375 fmult = min(FMAX,max(FMIN,fmult));
377 f[ntypes*ati+atj] *= fmult;
378 f[ntypes*atj+ati] *= fmult;
381 for(i=0; (i<ntypes); i++) {
382 f[ntypes*ati+i] *= fmult;
383 f[ntypes*i+ati] *= fmult;
390 static real calc_deviation(real xav,real xt,real x0)
392 /* This may prevent overshooting in GCT coupling... */
398 dev = min(xav-x0,xt-x0);
404 dev = max(xav-x0,xt-x0);
412 static real calc_dist(FILE *log,rvec x[])
414 static gmx_bool bFirst=TRUE;
415 static gmx_bool bDist;
421 if ((buf = getenv("DISTGCT")) == NULL)
424 bDist = (sscanf(buf,"%d%d",&i1,&i2) == 2);
426 fprintf(log,"GCT: Will couple to distance between %d and %d\n",i1,i2);
428 fprintf(log,"GCT: Will not couple to distances\n");
433 rvec_sub(x[i1],x[i2],dx);
440 real run_aver(real old,real cur,int step,int nmem)
444 return ((nmem-1)*old+cur)/nmem;
447 static void set_act_value(t_coupl_rec *tcr,int index,real val,int step)
449 tcr->act_value[index] = val;
450 tcr->av_value[index] = run_aver(tcr->av_value[index],val,step,tcr->nmemory);
453 static void upd_f_value(FILE *log,int atnr,real xi,real dt,real factor,
454 real ff[],int ati,int atj)
459 fff = 1 + (dt/xi) * factor;
461 set_factor_matrix(atnr,ff,sqrt(fff),ati,atj);
465 static void dump_fm(FILE *fp,int n,real f[],char *s)
469 fprintf(fp,"Factor matrix (all numbers -1) %s\n",s);
470 for(i=0; (i<n); i++) {
472 fprintf(fp," %10.3e",f[n*i+j]-1.0);
477 void do_coupling(FILE *log,const output_env_t oenv,int nfile,
478 const t_filenm fnm[], t_coupl_rec *tcr,real t,
479 int step,real ener[], t_forcerec *fr,t_inputrec *ir,
480 gmx_bool bMaster, t_mdatoms *md,t_idef *idef,real mu_aver,int nmols,
481 t_commrec *cr,matrix box,tensor virial,
482 tensor pres,rvec mu_tot,
483 rvec x[],rvec f[],gmx_bool bDoIt)
485 #define enm2Debye 48.0321
486 #define d2e(x) (x)/enm2Debye
487 #define enm2kjmol(x) (x)*0.0143952 /* = 2.0*4.0*M_PI*EPSILON0 */
489 static real *f6,*f12,*fa,*fb,*fc,*fq;
490 static gmx_bool bFirst = TRUE;
492 int i,j,ati,atj,atnr2,type,ftype;
493 real deviation[eoObsNR],prdev[eoObsNR],epot0,dist,rmsf;
494 real ff6,ff12,ffa,ffb,ffc,ffq,factor,dt,mu_ind;
495 real Epol,Eintern,Virial,muabs,xiH=-1,xiS=-1,xi6,xi12;
497 gmx_bool bTest,bPrint;
501 t_coupl_iparams *tip;
503 atnr2 = idef->atnr * idef->atnr;
506 fprintf(log,"GCT: this is parallel\n");
508 fprintf(log,"GCT: this is not parallel\n");
515 snew(fq, idef->atnr);
522 for(i=0; (i<ir->opts.ngtc); i++) {
523 nrdf += ir->opts.nrdf[i];
524 TTT += ir->opts.nrdf[i]*ir->opts.ref_t[i];
528 /* Calculate reference virial from reference temperature and pressure */
529 tcr->ref_value[eoVir] = 0.5*BOLTZ*nrdf*TTT - (3.0/2.0)*
530 Vol*tcr->ref_value[eoPres];
532 fprintf(log,"GCT: TTT = %g, nrdf = %d, vir0 = %g, Vol = %g\n",
533 TTT,nrdf,tcr->ref_value[eoVir],Vol);
539 bPrint = MASTER(cr) && do_per_step(step,ir->nstlog);
542 /* Initiate coupling to the reference pressure and temperature to start
546 for(i=0; (i<eoObsNR); i++)
547 tcr->av_value[i] = tcr->ref_value[i];
548 if ((tcr->ref_value[eoDipole]) != 0.0) {
549 mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
550 Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
551 tcr->av_value[eoEpot] -= Epol;
552 fprintf(log,"GCT: mu_aver = %g(D), mu_ind = %g(D), Epol = %g (kJ/mol)\n",
553 mu_aver*enm2Debye,mu_ind*enm2Debye,Epol);
557 /* We want to optimize the LJ params, usually to the Vaporization energy
558 * therefore we only count intermolecular degrees of freedom.
559 * Note that this is now optional. switch UseEinter to yes in your gct file
562 dist = calc_dist(log,x);
563 muabs = norm(mu_tot);
564 Eintern = Ecouple(tcr,ener)/nmols;
565 Virial = virial[XX][XX]+virial[YY][YY]+virial[ZZ][ZZ];
567 /*calc_force(md->nr,f,fmol);*/
570 /* Use a memory of tcr->nmemory steps, so we actually couple to the
571 * average observable over the last tcr->nmemory steps. This may help
572 * in avoiding local minima in parameter space.
574 set_act_value(tcr,eoPres, ener[F_PRES],step);
575 set_act_value(tcr,eoEpot, Eintern, step);
576 set_act_value(tcr,eoVir, Virial, step);
577 set_act_value(tcr,eoDist, dist, step);
578 set_act_value(tcr,eoMu, muabs, step);
579 set_act_value(tcr,eoFx, fmol[0][XX], step);
580 set_act_value(tcr,eoFy, fmol[0][YY], step);
581 set_act_value(tcr,eoFz, fmol[0][ZZ], step);
582 set_act_value(tcr,eoPx, pres[XX][XX],step);
583 set_act_value(tcr,eoPy, pres[YY][YY],step);
584 set_act_value(tcr,eoPz, pres[ZZ][ZZ],step);
586 epot0 = tcr->ref_value[eoEpot];
587 /* If dipole != 0.0 assume we want to use polarization corrected coupling */
588 if ((tcr->ref_value[eoDipole]) != 0.0) {
589 mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
591 Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
595 fprintf(debug,"mu_ind: %g (%g D) mu_aver: %g (%g D)\n",
596 mu_ind,mu_ind*enm2Debye,mu_aver,mu_aver*enm2Debye);
597 fprintf(debug,"Eref %g Epol %g Erunav %g Eact %g\n",
598 tcr->ref_value[eoEpot],Epol,tcr->av_value[eoEpot],
599 tcr->act_value[eoEpot]);
604 pr_ff(tcr,t,idef,cr,nfile,fnm,oenv);
606 /* Calculate the deviation of average value from the target value */
607 for(i=0; (i<eoObsNR); i++) {
608 deviation[i] = calc_deviation(tcr->av_value[i],tcr->act_value[i],
610 prdev[i] = tcr->ref_value[i] - tcr->act_value[i];
612 deviation[eoEpot] = calc_deviation(tcr->av_value[eoEpot],tcr->act_value[eoEpot],
614 prdev[eoEpot] = epot0 - tcr->act_value[eoEpot];
617 pr_dev(tcr,t,prdev,cr,nfile,fnm,oenv);
619 /* First set all factors to 1 */
620 for(i=0; (i<atnr2); i++) {
621 f6[i] = f12[i] = fa[i] = fb[i] = fc[i] = 1.0;
623 for(i=0; (i<idef->atnr); i++)
626 /* Now compute the actual coupling compononents */
629 for(i=0; (i<tcr->nLJ); i++) {
630 tclj=&(tcr->tcLJ[i]);
637 if (tclj->eObs == eoForce) {
638 gmx_fatal(FARGS,"Hack code for this to work again ");
640 fprintf(debug,"Have computed derivatives: xiH = %g, xiS = %g\n",xiH,xiS);
650 gmx_fatal(FARGS,"No H, no Shell, edit code at %s, line %d\n",
653 set_factor_matrix(idef->atnr,f6, sqrt(ff6), ati,atj);
655 set_factor_matrix(idef->atnr,f12,sqrt(ff12),ati,atj);
659 fprintf(debug,"tcr->tcLJ[%d].xi_6 = %g, xi_12 = %g deviation = %g\n",i,
660 tclj->xi_6,tclj->xi_12,deviation[tclj->eObs]);
661 factor=deviation[tclj->eObs];
663 upd_f_value(log,idef->atnr,tclj->xi_6, dt,factor,f6, ati,atj);
664 upd_f_value(log,idef->atnr,tclj->xi_12,dt,factor,f12,ati,atj);
672 dump_fm(log,idef->atnr,f6,"f6");
673 dump_fm(log,idef->atnr,f12,"f12");
676 upd_nbfplj(log,fr->nbfp,idef->atnr,f6,f12,tcr->combrule);
678 /* Copy for printing */
679 for(i=0; (i<tcr->nLJ); i++) {
680 tclj=&(tcr->tcLJ[i]);
687 /* nbfp now includes the 6.0/12.0 derivative prefactors */
688 tclj->c6 = C6(fr->nbfp,fr->ntype,ati,atj)/6.0;
689 tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj)/12.0;
694 for(i=0; (i<tcr->nBU); i++) {
695 tcbu = &(tcr->tcBU[i]);
696 factor = deviation[tcbu->eObs];
700 upd_f_value(log,idef->atnr,tcbu->xi_a,dt,factor,fa,ati,atj);
701 upd_f_value(log,idef->atnr,tcbu->xi_b,dt,factor,fb,ati,atj);
702 upd_f_value(log,idef->atnr,tcbu->xi_c,dt,factor,fc,ati,atj);
710 upd_nbfpbu(log,fr->nbfp,idef->atnr,fa,fb,fc);
711 /* Copy for printing */
712 for(i=0; (i<tcr->nBU); i++) {
713 tcbu=&(tcr->tcBU[i]);
720 /* nbfp now includes the 6.0 derivative prefactors */
721 tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
722 tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
723 tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj)/6.0;
725 fprintf(debug,"buck (type=%d) = %e, %e, %e\n",
726 tcbu->at_i,tcbu->a,tcbu->b,tcbu->c);
730 for(i=0; (i<tcr->nQ); i++) {
733 ffq = 1.0 + (dt/tcq->xi_Q) * deviation[tcq->eObs];
736 fq[tcq->at_i] *= ffq;
740 gprod(cr,idef->atnr,fq);
742 for(j=0; (j<md->nr); j++) {
743 md->chargeA[j] *= fq[md->typeA[j]];
745 for(i=0; (i<tcr->nQ); i++) {
747 for(j=0; (j<md->nr); j++) {
748 if (md->typeA[j] == tcq->at_i) {
749 tcq->Q = md->chargeA[j];
754 gmx_fatal(FARGS,"Coupling type %d not found",tcq->at_i);
756 for(i=0; (i<tcr->nIP); i++) {
757 tip = &(tcr->tIP[i]);
759 ftype = idef->functype[type];
760 factor = dt*deviation[tip->eObs];
764 if (tip->xi.harmonic.krA) idef->iparams[type].harmonic.krA *= (1+factor/tip->xi.harmonic.krA);
765 if (tip->xi.harmonic.rA) idef->iparams[type].harmonic.rA *= (1+factor/tip->xi.harmonic.rA);
770 tip->iprint=idef->iparams[type];