4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-2001
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
34 * Gromacs Runs On Most of All Computer Systems
36 static char *SRCID_gctio_c = "$Id$";
51 char *eoNames[eoNR] = {
52 "Pres", "Epot", "Vir", "Dist", "Mu", "Force", "Fx", "Fy", "Fz",
54 "Polarizability", "Dipole", "Memory", "UseEinter", "UseVirial"
57 static int Name2eo(char *s)
63 for(i=0; (i<eoNR); i++)
64 if (strcasecmp(s,eoNames[i]) == 0) {
66 fprintf(stderr,"Coupling to observable %d (%s)\n",res,eoNames[res]);
73 static char *NoYes[] = { "No", "Yes" };
75 static void send_tcr(int dest,t_coupl_rec *tcr)
77 nblocktx(dest,eoObsNR,tcr->ref_value);
78 blocktx(dest,tcr->nLJ);
79 nblocktx(dest,tcr->nLJ,tcr->tcLJ);
80 blocktx(dest,tcr->nBU);
81 nblocktx(dest,tcr->nBU,tcr->tcBU);
82 blocktx(dest,tcr->nQ);
83 nblocktx(dest,tcr->nQ,tcr->tcQ);
86 static void rec_tcr(int src,t_coupl_rec *tcr)
88 nblockrx(src,eoObsNR,tcr->ref_value);
90 blockrx(src,tcr->nLJ);
91 snew(tcr->tcLJ,tcr->nLJ);
92 nblockrx(src,tcr->nLJ,tcr->tcLJ);
94 blockrx(src,tcr->nBU);
95 snew(tcr->tcBU,tcr->nBU);
96 nblockrx(src,tcr->nBU,tcr->tcBU);
99 snew(tcr->tcQ,tcr->nQ);
100 nblockrx(src,tcr->nQ,tcr->tcQ);
103 void comm_tcr(FILE *log,t_commrec *cr,t_coupl_rec **tcr)
108 send_tcr(cr->left,*tcr);
109 rec_tcr(cr->right,&shit);
113 rec_tcr(cr->right,*tcr);
114 send_tcr(cr->left,*tcr);
118 static void clear_lj(t_coupl_LJ *tc)
130 static void clear_bu(t_coupl_BU *tc)
144 static void clear_q(t_coupl_Q *tc)
153 void copy_ff(t_coupl_rec *tcr,t_forcerec *fr,t_mdatoms *md,t_idef *idef)
155 int i,j,ati,atj,type;
160 /* Save values for printing */
161 for(i=0; (i<tcr->nLJ); i++) {
162 tclj = &(tcr->tcLJ[i]);
168 tclj->c6 = C6(fr->nbfp,fr->ntype,ati,atj);
169 tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj);
172 for(i=0; (i<tcr->nBU); i++) {
173 tcbu = &(tcr->tcBU[i]);
179 tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
180 tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
181 tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj);
184 for(i=0; (i<tcr->nQ); i++) {
185 tcq = &(tcr->tcQ[i]);
186 for(j=0; (j<md->nr); j++) {
187 if (md->typeA[j] == tcq->at_i) {
188 tcr->tcQ[i].Q = md->chargeA[j];
193 for(i=0; (i<tcr->nIP); i++) {
194 /* Let's just copy the whole struct !*/
195 type = tcr->tIP[i].type;
196 tcr->tIP[i].iprint=idef->iparams[type];
200 void write_gct(char *fn,t_coupl_rec *tcr,t_idef *idef)
207 fprintf(fp,"%-15s = %12g ; Reference pressure for coupling\n",
208 eoNames[eoPres],tcr->ref_value[eoPres]);
209 fprintf(fp,"%-15s = %12g ; Reference potential energy\n",
210 eoNames[eoEpot],tcr->ref_value[eoEpot]);
211 fprintf(fp,"%-15s = %12g ; Reference distance\n",
212 eoNames[eoDist],tcr->ref_value[eoDist]);
213 fprintf(fp,"%-15s = %12g ; Reference dipole\n",
214 eoNames[eoMu],tcr->ref_value[eoMu]);
215 fprintf(fp,"%-15s = %12g ; Reference force\n",
216 eoNames[eoForce],tcr->ref_value[eoForce]);
217 fprintf(fp,"%-15s = %12g ; Reference force in X dir\n",
218 eoNames[eoFx],tcr->ref_value[eoFx]);
219 fprintf(fp,"%-15s = %12g ; Reference force in Y dir\n",
220 eoNames[eoFy],tcr->ref_value[eoFy]);
221 fprintf(fp,"%-15s = %12g ; Reference force in Z dir\n",
222 eoNames[eoFz],tcr->ref_value[eoFz]);
223 fprintf(fp,"%-15s = %12g ; Reference pres in X dir\n",
224 eoNames[eoPx],tcr->ref_value[eoPx]);
225 fprintf(fp,"%-15s = %12g ; Reference pres in Y dir\n",
226 eoNames[eoPy],tcr->ref_value[eoPy]);
227 fprintf(fp,"%-15s = %12g ; Reference pres in Z dir\n",
228 eoNames[eoPz],tcr->ref_value[eoPz]);
229 fprintf(fp,"%-15s = %12g ; Polarizability used for the Epot correction\n",
230 eoNames[eoPolarizability],tcr->ref_value[eoPolarizability]);
231 fprintf(fp,"%-15s = %12g ; Gas phase dipole moment used for Epot correction\n",
232 eoNames[eoDipole],tcr->ref_value[eoDipole]);
233 fprintf(fp,"%-15s = %12d ; Memory for coupling. Makes it converge faster.\n",
234 eoNames[eoMemory],tcr->nmemory);
235 fprintf(fp,"%-15s = %12s ; Use intermolecular Epot only (LJ+Coul)\n",
236 eoNames[eoInter],yesno_names[tcr->bInter]);
237 fprintf(fp,"%-15s = %12s ; Use virial iso pressure\n",
238 eoNames[eoUseVirial],yesno_names[tcr->bVirial]);
240 fprintf(fp,"\n; Q-Coupling %6s %12s\n","type","xi");
241 for(i=0; (i<tcr->nQ); i++) {
242 fprintf(fp,"%-8s = %8s %6d %12g\n",
243 "Q",eoNames[tcr->tcQ[i].eObs],tcr->tcQ[i].at_i,tcr->tcQ[i].xi_Q);
246 fprintf(fp,"\n; %8s %8s %6s %6s %12s %12s\n","Couple","To",
247 "i-type","j-type","xi-c6","xi-c12");
248 fprintf(fp,"; j-type == -1 means mixing rules will be applied!\n");
249 for(i=0; (i<tcr->nLJ); i++) {
250 fprintf(fp,"%-8s = %8s %6d %6d %12g %12g\n",
251 "LJ",eoNames[tcr->tcLJ[i].eObs],
252 tcr->tcLJ[i].at_i,tcr->tcLJ[i].at_j,
253 tcr->tcLJ[i].xi_6,tcr->tcLJ[i].xi_12);
256 fprintf(fp,"\n; %8s %8s %6s %6s %12s %12s %12s\n","Couple","To",
257 "i-type","j-type","xi-A","xi-B","xi-C");
258 fprintf(fp,"; j-type == -1 means mixing rules will be applied!\n");
259 for(i=0; (i<tcr->nBU); i++) {
260 fprintf(fp,"%-8s = %8s %6d %6d %12g %12g %12g\n",
261 "BU",eoNames[tcr->tcBU[i].eObs],
262 tcr->tcBU[i].at_i,tcr->tcBU[i].at_j,
263 tcr->tcBU[i].xi_a,tcr->tcBU[i].xi_b,tcr->tcBU[i].xi_c);
266 fprintf(fp,"\n; More Coupling\n");
267 for(i=0; (i<tcr->nIP); i++) {
268 ftype=idef->functype[tcr->tIP[i].type];
271 fprintf(fp,"%-15s = %-8s %4d %12g %12g\n",
272 "Bonds",eoNames[tcr->tIP[i].eObs],tcr->tIP[i].type,
273 tcr->tIP[i].xi.harmonic.krA,
274 tcr->tIP[i].xi.harmonic.rA);
277 fprintf(stderr,"ftype %s not supported (yet)\n",
278 interaction_function[ftype].longname);
284 static bool add_lj(int *nLJ,t_coupl_LJ **tcLJ,char *s,bool bObsUsed[])
290 if (sscanf(s,"%s%d%d%lf%lf",buf,&ati,&atj,&xi6,&xi12) != 5)
292 if ((eo=Name2eo(buf)) == -1)
293 fatal_error(0,"Invalid observable for LJ coupling: %s",buf);
295 for(j=0; (j<*nLJ); j++) {
296 if ((((*tcLJ)[j].at_i == ati) && ((*tcLJ)[j].at_j == atj)) &&
297 ((*tcLJ)[j].xi_6 || (*tcLJ)[j].xi_12) &&
298 ((*tcLJ)[j].eObs == eo))
303 srenew((*tcLJ),*nLJ);
306 fprintf(stderr,"\n*** WARNING: overwriting entry for LJ coupling '%s'\n",s);
308 clear_lj(&((*tcLJ)[j]));
309 if (((*tcLJ)[j].eObs = eo) == -1) {
310 fatal_error(0,"Invalid observable for LJ coupling: %s",buf);
312 (*tcLJ)[j].at_i = ati;
313 (*tcLJ)[j].at_j = atj;
314 (*tcLJ)[j].xi_6 = xi6;
315 (*tcLJ)[j].xi_12 = xi12;
321 static bool add_bu(int *nBU,t_coupl_BU **tcBU,char *s,bool bObsUsed[])
327 if (sscanf(s,"%s%d%d%lf%lf%lf",buf,&ati,&atj,&xia,&xib,&xic) != 6)
329 if ((eo=Name2eo(buf)) == -1)
330 fatal_error(0,"Invalid observable for BU coupling: %s",buf);
332 for(j=0; (j<*nBU); j++) {
333 if ((((*tcBU)[j].at_i == ati) && ((*tcBU)[j].at_j == atj)) &&
334 ((*tcBU)[j].xi_a || (*tcBU)[j].xi_b || (*tcBU)[j].xi_c ) &&
335 ((*tcBU)[j].eObs == eo))
340 srenew((*tcBU),*nBU);
343 fprintf(stderr,"\n*** WARNING: overwriting entry for BU coupling '%s'\n",s);
345 clear_bu(&((*tcBU)[j]));
346 if (((*tcBU)[j].eObs = eo) == -1) {
347 fatal_error(0,"Invalid observable for BU coupling: %s",buf);
349 (*tcBU)[j].at_i = ati;
350 (*tcBU)[j].at_j = atj;
351 (*tcBU)[j].xi_a = xia;
352 (*tcBU)[j].xi_b = xib;
353 (*tcBU)[j].xi_c = xic;
359 static bool add_ip(int *nIP,t_coupl_iparams **tIP,char *s,int ftype,bool bObsUsed[])
367 /* Pick out the type */
368 if (sscanf(s,"%s%d",buf,&type) != 2)
370 if ((eo=Name2eo(buf)) == -1)
371 fatal_error(0,"Invalid observable for IP coupling: %s",buf);
373 /* Check whether this entry is there already */
374 for(i=0; (i<*nIP); i++) {
375 if ((*tIP)[i].type == type)
379 fprintf(stderr,"*** WARNING: overwriting entry for type %d\n",type);
386 if (sscanf(s,"%s%d%lf%lf",buf,&type,&kb,&b0) != 4)
390 (*tIP)[i].xi.harmonic.krA = kb;
391 (*tIP)[i].xi.harmonic.rA = b0;
395 fprintf(stderr,"ftype %s not supported (yet)\n",
396 interaction_function[ftype].longname);
402 static bool add_q(int *nQ,t_coupl_Q **tcQ,char *s,bool bObsUsed[])
408 if (sscanf(s,"%s%d%lf",buf,&ati,&xiQ) != 3)
411 for(j=0; (j<*nQ); j++) {
412 if ((*tcQ)[j].at_i == ati)
420 fprintf(stderr,"\n*** WARNING: overwriting entry for Q coupling '%s'\n",s);
422 clear_q(&((*tcQ)[j]));
423 eo = (*tcQ)[j].eObs = Name2eo(buf);
424 if ((*tcQ)[j].eObs == -1) {
425 fatal_error(0,"Invalid observable for Q coupling: %s",buf);
427 (*tcQ)[j].at_i = ati;
428 (*tcQ)[j].xi_Q = xiQ;
434 void read_gct(char *fn,t_coupl_rec *tcr)
437 int i,j,ninp,nQ,nLJ,nBU,nIP;
440 inp=read_inpfile(fn,&ninp);
441 for(i=0; (i<eoObsNR); i++) {
442 tcr->bObsUsed[i] = FALSE;
443 RTYPE (eoNames[i], tcr->ref_value[i], 0.0);
445 ITYPE (eoNames[eoMemory], tcr->nmemory, 1);
446 ETYPE (eoNames[eoInter], tcr->bInter, yesno_names);
447 ETYPE (eoNames[eoUseVirial], tcr->bVirial, yesno_names);
455 for(i=0; (i<ninp); i++) {
457 if (strcasecmp(inp[i].name,"LJ") == 0)
458 bWrong=add_lj(&nLJ,&(tcr->tcLJ),inp[i].value,tcr->bObsUsed);
459 else if (strcasecmp(inp[i].name,"BU") == 0)
460 bWrong=add_bu(&nBU,&(tcr->tcBU),inp[i].value,tcr->bObsUsed);
461 else if (strcasecmp(inp[i].name,"Q") == 0)
462 bWrong=add_q(&nQ,&(tcr->tcQ),inp[i].value,tcr->bObsUsed);
463 else if (strcasecmp(inp[i].name,"Bonds") == 0)
464 bWrong=add_ip(&nIP,&(tcr->tIP),inp[i].value,F_BONDS,tcr->bObsUsed);
467 fprintf(stderr,"Wrong line in %s: '%s = %s'\n",
468 fn,inp[i].name,inp[i].value);
469 /*sfree(inp[i].name);
470 sfree(inp[i].value);*/
472 /* Check which ones have to be printed */
473 for(i=1; (i<nQ); i++)
474 for(j=0; (j<i); j++) {
475 if (tcr->tcQ[i].at_i == tcr->tcQ[j].at_i)
476 tcr->tcQ[j].bPrint=FALSE;
478 for(i=1; (i<nLJ); i++)
479 for(j=0; (j<i); j++) {
480 if (((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_i) &&
481 (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_j)) ||
482 ((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_j) &&
483 (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_i)))
484 tcr->tcLJ[j].bPrint=FALSE;
487 for(i=1; (i<nBU); i++)
488 for(j=0; (j<i); j++) {
489 if (((tcr->tcBU[i].at_i == tcr->tcBU[j].at_i) &&
490 (tcr->tcBU[i].at_j == tcr->tcBU[j].at_j)) ||
491 ((tcr->tcBU[i].at_i == tcr->tcBU[j].at_j) &&
492 (tcr->tcBU[i].at_j == tcr->tcBU[j].at_i)))
493 tcr->tcBU[j].bPrint=FALSE;