3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
53 const char *eoNames[eoNR] = {
54 "Pres", "Epot", "Vir", "Dist", "Mu", "Force", "Fx", "Fy", "Fz",
56 "Polarizability", "Dipole", "Memory", "UseEinter", "UseVirial",
60 static int Name2eo(char *s)
66 for(i=0; (i<eoNR); i++)
67 if (gmx_strcasecmp(s,eoNames[i]) == 0) {
69 fprintf(stderr,"Coupling to observable %d (%s)\n",res,eoNames[res]);
76 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
77 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
78 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
80 static void low_comm_tcr(t_commrec *cr,t_coupl_rec *tcr)
82 nblock_bc(cr,eoObsNR,tcr->ref_value);
84 block_bc(cr,tcr->nLJ);
85 snew_bc(cr,tcr->tcLJ,tcr->nLJ);
86 nblock_bc(cr,tcr->nLJ,tcr->tcLJ);
88 block_bc(cr,tcr->nBU);
89 snew_bc(cr,tcr->tcBU,tcr->nBU);
90 nblock_bc(cr,tcr->nBU,tcr->tcBU);
93 snew_bc(cr,tcr->tcQ,tcr->nQ);
94 nblock_bc(cr,tcr->nQ,tcr->tcQ);
96 block_bc(cr,tcr->nmemory);
97 block_bc(cr,tcr->bInter);
98 block_bc(cr,tcr->bVirial);
99 block_bc(cr,tcr->combrule);
102 void comm_tcr(FILE *log,t_commrec *cr,t_coupl_rec **tcr)
107 low_comm_tcr(cr,*tcr);
110 static void clear_lj(t_coupl_LJ *tc)
122 static void clear_bu(t_coupl_BU *tc)
136 static void clear_q(t_coupl_Q *tc)
145 void copy_ff(t_coupl_rec *tcr,t_forcerec *fr,t_mdatoms *md,t_idef *idef)
147 int i,j,ati,atj,type;
152 /* Save values for printing */
153 for(i=0; (i<tcr->nLJ); i++) {
154 tclj = &(tcr->tcLJ[i]);
160 tclj->c6 = C6(fr->nbfp,fr->ntype,ati,atj);
161 tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj);
164 for(i=0; (i<tcr->nBU); i++) {
165 tcbu = &(tcr->tcBU[i]);
171 tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
172 tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
173 tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj);
176 for(i=0; (i<tcr->nQ); i++) {
177 tcq = &(tcr->tcQ[i]);
178 for(j=0; (j<md->nr); j++) {
179 if (md->typeA[j] == tcq->at_i) {
180 tcr->tcQ[i].Q = md->chargeA[j];
185 for(i=0; (i<tcr->nIP); i++) {
186 /* Let's just copy the whole struct !*/
187 type = tcr->tIP[i].type;
188 tcr->tIP[i].iprint=idef->iparams[type];
192 void write_gct(const char *fn,t_coupl_rec *tcr,t_idef *idef)
197 fp=gmx_fio_fopen(fn,"w");
199 fprintf(fp,"%-15s = %12g ; Reference pressure for coupling\n",
200 eoNames[eoPres],tcr->ref_value[eoPres]);
201 fprintf(fp,"%-15s = %12g ; Reference potential energy\n",
202 eoNames[eoEpot],tcr->ref_value[eoEpot]);
203 fprintf(fp,"%-15s = %12g ; Reference distance\n",
204 eoNames[eoDist],tcr->ref_value[eoDist]);
205 fprintf(fp,"%-15s = %12g ; Reference dipole\n",
206 eoNames[eoMu],tcr->ref_value[eoMu]);
207 fprintf(fp,"%-15s = %12g ; Reference force\n",
208 eoNames[eoForce],tcr->ref_value[eoForce]);
209 fprintf(fp,"%-15s = %12g ; Reference force in X dir\n",
210 eoNames[eoFx],tcr->ref_value[eoFx]);
211 fprintf(fp,"%-15s = %12g ; Reference force in Y dir\n",
212 eoNames[eoFy],tcr->ref_value[eoFy]);
213 fprintf(fp,"%-15s = %12g ; Reference force in Z dir\n",
214 eoNames[eoFz],tcr->ref_value[eoFz]);
215 fprintf(fp,"%-15s = %12g ; Reference pres in X dir\n",
216 eoNames[eoPx],tcr->ref_value[eoPx]);
217 fprintf(fp,"%-15s = %12g ; Reference pres in Y dir\n",
218 eoNames[eoPy],tcr->ref_value[eoPy]);
219 fprintf(fp,"%-15s = %12g ; Reference pres in Z dir\n",
220 eoNames[eoPz],tcr->ref_value[eoPz]);
221 fprintf(fp,"%-15s = %12g ; Polarizability used for the Epot correction\n",
222 eoNames[eoPolarizability],tcr->ref_value[eoPolarizability]);
223 fprintf(fp,"%-15s = %12g ; Gas phase dipole moment used for Epot correction\n",
224 eoNames[eoDipole],tcr->ref_value[eoDipole]);
225 fprintf(fp,"%-15s = %12d ; Memory for coupling. Makes it converge faster.\n",
226 eoNames[eoMemory],tcr->nmemory);
227 fprintf(fp,"%-15s = %12s ; Use intermolecular Epot only (LJ+Coul)\n",
228 eoNames[eoInter],yesno_names[tcr->bInter]);
229 fprintf(fp,"%-15s = %12s ; Use virial iso pressure\n",
230 eoNames[eoUseVirial],yesno_names[tcr->bVirial]);
231 fprintf(fp,"%-15s = %12d ; Combination rule, same coding as in grompp.\n",
232 eoNames[eoCombRule],tcr->combrule);
234 fprintf(fp,"\n; Q-Coupling %6s %12s\n","type","xi");
235 for(i=0; (i<tcr->nQ); i++) {
236 fprintf(fp,"%-8s = %8s %6d %12g\n",
237 "Q",eoNames[tcr->tcQ[i].eObs],tcr->tcQ[i].at_i,tcr->tcQ[i].xi_Q);
240 fprintf(fp,"\n; %8s %8s %6s %6s %12s %12s\n","Couple","To",
241 "i-type","j-type","xi-c6","xi-c12");
242 fprintf(fp,"; j-type == -1 means mixing rules will be applied!\n");
243 for(i=0; (i<tcr->nLJ); i++) {
244 fprintf(fp,"%-8s = %8s %6d %6d %12g %12g\n",
245 "LJ",eoNames[tcr->tcLJ[i].eObs],
246 tcr->tcLJ[i].at_i,tcr->tcLJ[i].at_j,
247 tcr->tcLJ[i].xi_6,tcr->tcLJ[i].xi_12);
250 fprintf(fp,"\n; %8s %8s %6s %6s %12s %12s %12s\n","Couple","To",
251 "i-type","j-type","xi-A","xi-B","xi-C");
252 fprintf(fp,"; j-type == -1 means mixing rules will be applied!\n");
253 for(i=0; (i<tcr->nBU); i++) {
254 fprintf(fp,"%-8s = %8s %6d %6d %12g %12g %12g\n",
255 "BU",eoNames[tcr->tcBU[i].eObs],
256 tcr->tcBU[i].at_i,tcr->tcBU[i].at_j,
257 tcr->tcBU[i].xi_a,tcr->tcBU[i].xi_b,tcr->tcBU[i].xi_c);
260 fprintf(fp,"\n; More Coupling\n");
261 for(i=0; (i<tcr->nIP); i++) {
262 ftype=idef->functype[tcr->tIP[i].type];
265 fprintf(fp,"%-15s = %-8s %4d %12g %12g\n",
266 "Bonds",eoNames[tcr->tIP[i].eObs],tcr->tIP[i].type,
267 tcr->tIP[i].xi.harmonic.krA,
268 tcr->tIP[i].xi.harmonic.rA);
271 fprintf(stderr,"ftype %s not supported (yet)\n",
272 interaction_function[ftype].longname);
278 static gmx_bool add_lj(int *nLJ,t_coupl_LJ **tcLJ,char *s,gmx_bool bObsUsed[])
284 if (sscanf(s,"%s%d%d%lf%lf",buf,&ati,&atj,&xi6,&xi12) != 5)
286 if ((eo=Name2eo(buf)) == -1)
287 gmx_fatal(FARGS,"Invalid observable for LJ coupling: %s",buf);
289 for(j=0; (j<*nLJ); j++) {
290 if ((((*tcLJ)[j].at_i == ati) && ((*tcLJ)[j].at_j == atj)) &&
291 ((*tcLJ)[j].xi_6 || (*tcLJ)[j].xi_12) &&
292 ((*tcLJ)[j].eObs == eo))
297 srenew((*tcLJ),*nLJ);
300 fprintf(stderr,"\n*** WARNING: overwriting entry for LJ coupling '%s'\n",s);
302 clear_lj(&((*tcLJ)[j]));
303 if (((*tcLJ)[j].eObs = eo) == -1) {
304 gmx_fatal(FARGS,"Invalid observable for LJ coupling: %s",buf);
306 (*tcLJ)[j].at_i = ati;
307 (*tcLJ)[j].at_j = atj;
308 (*tcLJ)[j].xi_6 = xi6;
309 (*tcLJ)[j].xi_12 = xi12;
315 static gmx_bool add_bu(int *nBU,t_coupl_BU **tcBU,char *s,gmx_bool bObsUsed[])
321 if (sscanf(s,"%s%d%d%lf%lf%lf",buf,&ati,&atj,&xia,&xib,&xic) != 6)
323 if ((eo=Name2eo(buf)) == -1)
324 gmx_fatal(FARGS,"Invalid observable for BU coupling: %s",buf);
326 for(j=0; (j<*nBU); j++) {
327 if ((((*tcBU)[j].at_i == ati) && ((*tcBU)[j].at_j == atj)) &&
328 ((*tcBU)[j].xi_a || (*tcBU)[j].xi_b || (*tcBU)[j].xi_c ) &&
329 ((*tcBU)[j].eObs == eo))
334 srenew((*tcBU),*nBU);
337 fprintf(stderr,"\n*** WARNING: overwriting entry for BU coupling '%s'\n",s);
339 clear_bu(&((*tcBU)[j]));
340 if (((*tcBU)[j].eObs = eo) == -1) {
341 gmx_fatal(FARGS,"Invalid observable for BU coupling: %s",buf);
343 (*tcBU)[j].at_i = ati;
344 (*tcBU)[j].at_j = atj;
345 (*tcBU)[j].xi_a = xia;
346 (*tcBU)[j].xi_b = xib;
347 (*tcBU)[j].xi_c = xic;
353 static gmx_bool add_ip(int *nIP,t_coupl_iparams **tIP,char *s,int ftype,gmx_bool bObsUsed[])
361 /* Pick out the type */
362 if (sscanf(s,"%s%d",buf,&type) != 2)
364 if ((eo=Name2eo(buf)) == -1)
365 gmx_fatal(FARGS,"Invalid observable for IP coupling: %s",buf);
367 /* Check whether this entry is there already */
368 for(i=0; (i<*nIP); i++) {
369 if ((*tIP)[i].type == type)
373 fprintf(stderr,"*** WARNING: overwriting entry for type %d\n",type);
380 if (sscanf(s,"%s%d%lf%lf",buf,&type,&kb,&b0) != 4)
384 (*tIP)[i].xi.harmonic.krA = kb;
385 (*tIP)[i].xi.harmonic.rA = b0;
389 fprintf(stderr,"ftype %s not supported (yet)\n",
390 interaction_function[ftype].longname);
396 static gmx_bool add_q(int *nQ,t_coupl_Q **tcQ,char *s,gmx_bool bObsUsed[])
402 if (sscanf(s,"%s%d%lf",buf,&ati,&xiQ) != 3)
405 for(j=0; (j<*nQ); j++) {
406 if ((*tcQ)[j].at_i == ati)
414 fprintf(stderr,"\n*** WARNING: overwriting entry for Q coupling '%s'\n",s);
416 clear_q(&((*tcQ)[j]));
417 eo = (*tcQ)[j].eObs = Name2eo(buf);
418 if ((*tcQ)[j].eObs == -1) {
419 gmx_fatal(FARGS,"Invalid observable for Q coupling: %s",buf);
421 (*tcQ)[j].at_i = ati;
422 (*tcQ)[j].xi_Q = xiQ;
428 void read_gct(const char *fn,t_coupl_rec *tcr)
432 int i,j,ninp,nQ,nLJ,nBU,nIP;
435 wi = init_warning(FALSE,0);
437 inp=read_inpfile(fn,&ninp,NULL,wi);
439 for(i=0; (i<eoObsNR); i++) {
440 tcr->bObsUsed[i] = FALSE;
441 RTYPE (eoNames[i], tcr->ref_value[i], 0.0);
443 ITYPE (eoNames[eoMemory], tcr->nmemory, 1);
444 ETYPE (eoNames[eoInter], tcr->bInter, yesno_names);
445 ETYPE (eoNames[eoUseVirial], tcr->bVirial, yesno_names);
446 ITYPE (eoNames[eoCombRule], tcr->combrule, 1);
453 for(i=0; (i<ninp); i++) {
455 if (gmx_strcasecmp(inp[i].name,"LJ") == 0)
456 bWrong=add_lj(&nLJ,&(tcr->tcLJ),inp[i].value,tcr->bObsUsed);
457 else if (gmx_strcasecmp(inp[i].name,"BU") == 0)
458 bWrong=add_bu(&nBU,&(tcr->tcBU),inp[i].value,tcr->bObsUsed);
459 else if (gmx_strcasecmp(inp[i].name,"Q") == 0)
460 bWrong=add_q(&nQ,&(tcr->tcQ),inp[i].value,tcr->bObsUsed);
461 else if (gmx_strcasecmp(inp[i].name,"Bonds") == 0)
462 bWrong=add_ip(&nIP,&(tcr->tIP),inp[i].value,F_BONDS,tcr->bObsUsed);
465 fprintf(stderr,"Wrong line in %s: '%s = %s'\n",
466 fn,inp[i].name,inp[i].value);
467 /*sfree(inp[i].name);
468 sfree(inp[i].value);*/
470 /* Check which ones have to be printed */
471 for(i=1; (i<nQ); i++)
472 for(j=0; (j<i); j++) {
473 if (tcr->tcQ[i].at_i == tcr->tcQ[j].at_i)
474 tcr->tcQ[j].bPrint=FALSE;
476 for(i=1; (i<nLJ); i++)
477 for(j=0; (j<i); j++) {
478 if (((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_i) &&
479 (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_j)) ||
480 ((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_j) &&
481 (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_i)))
482 tcr->tcLJ[j].bPrint=FALSE;
485 for(i=1; (i<nBU); i++)
486 for(j=0; (j<i); j++) {
487 if (((tcr->tcBU[i].at_i == tcr->tcBU[j].at_i) &&
488 (tcr->tcBU[i].at_j == tcr->tcBU[j].at_j)) ||
489 ((tcr->tcBU[i].at_i == tcr->tcBU[j].at_j) &&
490 (tcr->tcBU[i].at_j == tcr->tcBU[j].at_i)))
491 tcr->tcBU[j].bPrint=FALSE;
501 done_warning(wi,FARGS);