src/kernel/gctio.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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,2013, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "typedefs.h"
  43 #include "xmdrun.h"
  44 #include "futil.h"
  45 #include "xvgr.h"
  46 #include "macros.h"
  47 #include "physics.h"
  48 #include "network.h"
  49 #include "smalloc.h"
  50 #include "string2.h"
  51 #include "readinp.h"
  52 #include "readir.h"
  53 #include "filenm.h"
  54 #include "names.h"
  55 #include "gmxfio.h"
  56
  57 const char *eoNames[eoNR] = {
  58   "Pres", "Epot", "Vir", "Dist", "Mu", "Force", "Fx", "Fy", "Fz",
  59   "Px", "Py", "Pz",
  60   "Polarizability", "Dipole", "Memory", "UseEinter", "UseVirial",
  61   "CombinationRule"
  62 };
  63
  64 static int Name2eo(char *s)
  65 {
  66   int i,res;
  67
  68   res=-1;
  69
  70   for(i=0; (i<eoNR); i++)
  71     if (gmx_strcasecmp(s,eoNames[i]) == 0) {
  72       res=i;
  73       fprintf(stderr,"Coupling to observable %d (%s)\n",res,eoNames[res]);
  74       break;
  75     }
  76
  77   return res;
  78 }
  79
  80 #define  block_bc(cr,   d) gmx_bcast(     sizeof(d),     &(d),(cr))
  81 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
  82 #define   snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
  83
  84 static void low_comm_tcr(t_commrec *cr,t_coupl_rec *tcr)
  85 {
  86   nblock_bc(cr,eoObsNR,tcr->ref_value);
  87
  88   block_bc(cr,tcr->nLJ);
  89   snew_bc(cr,tcr->tcLJ,tcr->nLJ);
  90   nblock_bc(cr,tcr->nLJ,tcr->tcLJ);
  91
  92   block_bc(cr,tcr->nBU);
  93   snew_bc(cr,tcr->tcBU,tcr->nBU);
  94   nblock_bc(cr,tcr->nBU,tcr->tcBU);
  95
  96   block_bc(cr,tcr->nQ);
  97   snew_bc(cr,tcr->tcQ,tcr->nQ);
  98   nblock_bc(cr,tcr->nQ,tcr->tcQ);
  99
 100   block_bc(cr,tcr->nmemory);
 101   block_bc(cr,tcr->bInter);
 102   block_bc(cr,tcr->bVirial);
 103   block_bc(cr,tcr->combrule);
 104 }
 105
 106 void comm_tcr(FILE *log,t_commrec *cr,t_coupl_rec **tcr)
 107 {
 108   if (!MASTER(cr))
 109     snew(*tcr,1);
 110
 111   low_comm_tcr(cr,*tcr);
 112 }
 113
 114 static void clear_lj(t_coupl_LJ *tc)
 115 {
 116   tc->at_i   = 0;
 117   tc->at_j   = 0;
 118   tc->eObs   = -1;
 119   tc->bPrint = TRUE;
 120   tc->c6     = 0.0;
 121   tc->c12    = 0.0;
 122   tc->xi_6   = 0.0;
 123   tc->xi_12  = 0.0;
 124 }
 125
 126 static void clear_bu(t_coupl_BU *tc)
 127 {
 128   tc->at_i   = 0;
 129   tc->at_j   = 0;
 130   tc->eObs   = -1;
 131   tc->bPrint = TRUE;
 132   tc->a      = 0.0;
 133   tc->b      = 0.0;
 134   tc->c      = 0.0;
 135   tc->xi_a   = 0.0;
 136   tc->xi_b   = 0.0;
 137   tc->xi_c   = 0.0;
 138 }
 139
 140 static void clear_q(t_coupl_Q *tc)
 141 {
 142   tc->at_i   = 0;
 143   tc->eObs   = -1;
 144   tc->bPrint = TRUE;
 145   tc->Q      = 0.0;
 146   tc->xi_Q   = 0.0;
 147 }
 148
 149 void copy_ff(t_coupl_rec *tcr,t_forcerec *fr,t_mdatoms *md,t_idef *idef)
 150 {
 151   int        i,j,ati,atj,type;
 152   t_coupl_LJ *tclj;
 153   t_coupl_BU *tcbu;
 154   t_coupl_Q  *tcq;
 155
 156   /* Save values for printing */
 157   for(i=0; (i<tcr->nLJ); i++) {
 158     tclj = &(tcr->tcLJ[i]);
 159
 160     ati  = tclj->at_i;
 161     atj  = tclj->at_j;
 162     if (atj == -1)
 163       atj = ati;
 164     /* nbfp now includes the 6.0/12.0 derivative prefactors */
 165     tclj->c6  = C6(fr->nbfp,fr->ntype,ati,atj)/6.0;
 166     tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj)/12.0;
 167   }
 168
 169   for(i=0; (i<tcr->nBU); i++) {
 170     tcbu = &(tcr->tcBU[i]);
 171
 172     ati  = tcbu->at_i;
 173     atj  = tcbu->at_j;
 174     if (atj == -1)
 175       atj = ati;
 176     /* nbfp now includes the 6.0 derivative prefactor */
 177     tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
 178     tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
 179     tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj)/6.0;
 180   }
 181
 182   for(i=0; (i<tcr->nQ); i++) {
 183     tcq = &(tcr->tcQ[i]);
 184     for(j=0; (j<md->nr); j++) {
 185       if (md->typeA[j] == tcq->at_i) {
 186         tcr->tcQ[i].Q = md->chargeA[j];
 187         break;
 188       }
 189     }
 190   }
 191   for(i=0; (i<tcr->nIP); i++) {
 192     /* Let's just copy the whole struct !*/
 193     type = tcr->tIP[i].type;
 194     tcr->tIP[i].iprint=idef->iparams[type];
 195   }
 196 }
 197
 198 void write_gct(const char *fn,t_coupl_rec *tcr,t_idef *idef)
 199 {
 200   FILE *fp;
 201   int  i,ftype;
 202
 203   fp=gmx_fio_fopen(fn,"w");
 204   nice_header(fp,fn);
 205   fprintf(fp,"%-15s = %12g  ; Reference pressure for coupling\n",
 206           eoNames[eoPres],tcr->ref_value[eoPres]);
 207   fprintf(fp,"%-15s = %12g  ; Reference potential energy\n",
 208           eoNames[eoEpot],tcr->ref_value[eoEpot]);
 209   fprintf(fp,"%-15s = %12g  ; Reference distance\n",
 210           eoNames[eoDist],tcr->ref_value[eoDist]);
 211   fprintf(fp,"%-15s = %12g  ; Reference dipole\n",
 212           eoNames[eoMu],tcr->ref_value[eoMu]);
 213   fprintf(fp,"%-15s = %12g  ; Reference force\n",
 214           eoNames[eoForce],tcr->ref_value[eoForce]);
 215   fprintf(fp,"%-15s = %12g  ; Reference force in X dir\n",
 216           eoNames[eoFx],tcr->ref_value[eoFx]);
 217   fprintf(fp,"%-15s = %12g  ; Reference force in Y dir\n",
 218           eoNames[eoFy],tcr->ref_value[eoFy]);
 219   fprintf(fp,"%-15s = %12g  ; Reference force in Z dir\n",
 220           eoNames[eoFz],tcr->ref_value[eoFz]);
 221   fprintf(fp,"%-15s = %12g  ; Reference pres in X dir\n",
 222           eoNames[eoPx],tcr->ref_value[eoPx]);
 223   fprintf(fp,"%-15s = %12g  ; Reference pres in Y dir\n",
 224           eoNames[eoPy],tcr->ref_value[eoPy]);
 225   fprintf(fp,"%-15s = %12g  ; Reference pres in Z dir\n",
 226           eoNames[eoPz],tcr->ref_value[eoPz]);
 227   fprintf(fp,"%-15s = %12g  ; Polarizability used for the Epot correction\n",
 228           eoNames[eoPolarizability],tcr->ref_value[eoPolarizability]);
 229   fprintf(fp,"%-15s = %12g  ; Gas phase dipole moment used for Epot correction\n",
 230           eoNames[eoDipole],tcr->ref_value[eoDipole]);
 231   fprintf(fp,"%-15s = %12d  ; Memory for coupling. Makes it converge faster.\n",
 232           eoNames[eoMemory],tcr->nmemory);
 233   fprintf(fp,"%-15s = %12s  ; Use intermolecular Epot only (LJ+Coul)\n",
 234           eoNames[eoInter],yesno_names[tcr->bInter]);
 235   fprintf(fp,"%-15s = %12s  ; Use virial iso pressure\n",
 236           eoNames[eoUseVirial],yesno_names[tcr->bVirial]);
 237   fprintf(fp,"%-15s = %12d  ; Combination rule, same coding as in grompp.\n",
 238           eoNames[eoCombRule],tcr->combrule);
 239
 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);
 244   }
 245
 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);
 254   }
 255
 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);
 264   }
 265
 266   fprintf(fp,"\n; More Coupling\n");
 267   for(i=0; (i<tcr->nIP); i++) {
 268     ftype=idef->functype[tcr->tIP[i].type];
 269     switch (ftype) {
 270     case F_BONDS:
 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);
 275       break;
 276     default:
 277       fprintf(stderr,"ftype %s not supported (yet)\n",
 278               interaction_function[ftype].longname);
 279     }
 280   }
 281   gmx_fio_fclose(fp);
 282 }
 283
 284 static gmx_bool add_lj(int *nLJ,t_coupl_LJ **tcLJ,char *s,gmx_bool bObsUsed[])
 285 {
 286   int       j,ati,atj,eo;
 287   char      buf[256];
 288   double    xi6,xi12;
 289
 290   if (sscanf(s,"%s%d%d%lf%lf",buf,&ati,&atj,&xi6,&xi12) != 5)
 291     return TRUE;
 292   if ((eo=Name2eo(buf)) == -1)
 293     gmx_fatal(FARGS,"Invalid observable for LJ coupling: %s",buf);
 294
 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))
 299       break;
 300   }
 301   if (j == *nLJ) {
 302     ++(*nLJ);
 303     srenew((*tcLJ),*nLJ);
 304   }
 305   else
 306     fprintf(stderr,"\n*** WARNING: overwriting entry for LJ coupling '%s'\n",s);
 307
 308   clear_lj(&((*tcLJ)[j]));
 309   if (((*tcLJ)[j].eObs = eo) == -1) {
 310     gmx_fatal(FARGS,"Invalid observable for LJ coupling: %s",buf);
 311   }
 312   (*tcLJ)[j].at_i   = ati;
 313   (*tcLJ)[j].at_j   = atj;
 314   (*tcLJ)[j].xi_6   = xi6;
 315   (*tcLJ)[j].xi_12  = xi12;
 316   bObsUsed[eo] = TRUE;
 317
 318   return FALSE;
 319 }
 320
 321 static gmx_bool add_bu(int *nBU,t_coupl_BU **tcBU,char *s,gmx_bool bObsUsed[])
 322 {
 323   int       j,ati,atj,eo;
 324   char      buf[256];
 325   double    xia,xib,xic;
 326
 327   if (sscanf(s,"%s%d%d%lf%lf%lf",buf,&ati,&atj,&xia,&xib,&xic) != 6)
 328     return TRUE;
 329   if ((eo=Name2eo(buf)) == -1)
 330     gmx_fatal(FARGS,"Invalid observable for BU coupling: %s",buf);
 331
 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))
 336       break;
 337   }
 338   if (j == *nBU) {
 339     ++(*nBU);
 340     srenew((*tcBU),*nBU);
 341   }
 342   else
 343     fprintf(stderr,"\n*** WARNING: overwriting entry for BU coupling '%s'\n",s);
 344
 345   clear_bu(&((*tcBU)[j]));
 346   if (((*tcBU)[j].eObs = eo) == -1) {
 347     gmx_fatal(FARGS,"Invalid observable for BU coupling: %s",buf);
 348   }
 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;
 354   bObsUsed[eo] = TRUE;
 355
 356   return FALSE;
 357 }
 358
 359 static gmx_bool add_ip(int *nIP,t_coupl_iparams **tIP,char *s,int ftype,gmx_bool bObsUsed[])
 360 {
 361   int    i,eo,type;
 362   char   buf[256];
 363   double kb,b0;
 364
 365   switch (ftype) {
 366   case F_BONDS:
 367     /* Pick out the type */
 368     if (sscanf(s,"%s%d",buf,&type) != 2)
 369       return TRUE;
 370     if ((eo=Name2eo(buf)) == -1)
 371       gmx_fatal(FARGS,"Invalid observable for IP coupling: %s",buf);
 372
 373     /* Check whether this entry is there already */
 374     for(i=0; (i<*nIP); i++) {
 375       if ((*tIP)[i].type == type)
 376         break;
 377     }
 378     if (i < *nIP) {
 379       fprintf(stderr,"*** WARNING: overwriting entry for type %d\n",type);
 380     }
 381     else {
 382       i=*nIP;
 383       srenew((*tIP),i+1);
 384       (*nIP)++;
 385     }
 386     if (sscanf(s,"%s%d%lf%lf",buf,&type,&kb,&b0) != 4)
 387       return TRUE;
 388     (*tIP)[i].type=type;
 389     (*tIP)[i].eObs=eo;
 390     (*tIP)[i].xi.harmonic.krA = kb;
 391     (*tIP)[i].xi.harmonic.rA  = b0;
 392     bObsUsed[eo] = TRUE;
 393     break;
 394   default:
 395     fprintf(stderr,"ftype %s not supported (yet)\n",
 396             interaction_function[ftype].longname);
 397     return TRUE;
 398   }
 399   return FALSE;
 400 }
 401
 402 static gmx_bool add_q(int *nQ,t_coupl_Q **tcQ,char *s,gmx_bool bObsUsed[])
 403 {
 404   int       j,ati,eo;
 405   char      buf[256];
 406   double    xiQ;
 407
 408   if (sscanf(s,"%s%d%lf",buf,&ati,&xiQ) != 3)
 409     return TRUE;
 410
 411   for(j=0; (j<*nQ); j++) {
 412     if ((*tcQ)[j].at_i == ati)
 413       break;
 414   }
 415   if (j == *nQ) {
 416     ++(*nQ);
 417     srenew((*tcQ),*nQ);
 418   }
 419   else
 420     fprintf(stderr,"\n*** WARNING: overwriting entry for Q coupling '%s'\n",s);
 421
 422   clear_q(&((*tcQ)[j]));
 423   eo = (*tcQ)[j].eObs = Name2eo(buf);
 424   if ((*tcQ)[j].eObs == -1) {
 425     gmx_fatal(FARGS,"Invalid observable for Q coupling: %s",buf);
 426   }
 427   (*tcQ)[j].at_i   = ati;
 428   (*tcQ)[j].xi_Q  = xiQ;
 429   bObsUsed[eo] = TRUE;
 430
 431   return FALSE;
 432 }
 433
 434 void read_gct(const char *fn,t_coupl_rec *tcr)
 435 {
 436   warninp_t wi;
 437   t_inpfile *inp;
 438   int       i,j,ninp,nQ,nLJ,nBU,nIP;
 439   gmx_bool      bWrong;
 440
 441   wi = init_warning(FALSE,0);
 442
 443   inp=read_inpfile(fn,&ninp,NULL,wi);
 444
 445   for(i=0; (i<eoObsNR); i++) {
 446     tcr->bObsUsed[i] = FALSE;
 447     RTYPE (eoNames[i],  tcr->ref_value[i],      0.0);
 448   }
 449   ITYPE (eoNames[eoMemory],     tcr->nmemory,   1);
 450   ETYPE (eoNames[eoInter],      tcr->bInter,    yesno_names);
 451   ETYPE (eoNames[eoUseVirial],  tcr->bVirial,   yesno_names);
 452   ITYPE (eoNames[eoCombRule],   tcr->combrule,  1);
 453   tcr->tcLJ=NULL;
 454   tcr->tcBU=NULL;
 455   tcr->tcQ=NULL;
 456   tcr->tIP=NULL;
 457   nQ=nLJ=nBU=nIP=0;
 458
 459   for(i=0; (i<ninp); i++) {
 460     bWrong=FALSE;
 461     if (gmx_strcasecmp(inp[i].name,"LJ") == 0)
 462       bWrong=add_lj(&nLJ,&(tcr->tcLJ),inp[i].value,tcr->bObsUsed);
 463     else if (gmx_strcasecmp(inp[i].name,"BU") == 0)
 464       bWrong=add_bu(&nBU,&(tcr->tcBU),inp[i].value,tcr->bObsUsed);
 465     else if (gmx_strcasecmp(inp[i].name,"Q") == 0)
 466       bWrong=add_q(&nQ,&(tcr->tcQ),inp[i].value,tcr->bObsUsed);
 467     else if (gmx_strcasecmp(inp[i].name,"Bonds") == 0)
 468       bWrong=add_ip(&nIP,&(tcr->tIP),inp[i].value,F_BONDS,tcr->bObsUsed);
 469
 470     if (bWrong)
 471       fprintf(stderr,"Wrong line in %s: '%s = %s'\n",
 472               fn,inp[i].name,inp[i].value);
 473     /*sfree(inp[i].name);
 474       sfree(inp[i].value);*/
 475   }
 476   /* Check which ones have to be printed */
 477   for(i=1; (i<nQ); i++)
 478     for(j=0; (j<i); j++) {
 479       if (tcr->tcQ[i].at_i == tcr->tcQ[j].at_i)
 480         tcr->tcQ[j].bPrint=FALSE;
 481     }
 482   for(i=1; (i<nLJ); i++)
 483     for(j=0; (j<i); j++) {
 484       if (((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_i) &&
 485            (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_j)) ||
 486           ((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_j) &&
 487            (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_i)))
 488         tcr->tcLJ[j].bPrint=FALSE;
 489     }
 490
 491   for(i=1; (i<nBU); i++)
 492     for(j=0; (j<i); j++) {
 493       if (((tcr->tcBU[i].at_i == tcr->tcBU[j].at_i) &&
 494            (tcr->tcBU[i].at_j == tcr->tcBU[j].at_j)) ||
 495           ((tcr->tcBU[i].at_i == tcr->tcBU[j].at_j) &&
 496            (tcr->tcBU[i].at_j == tcr->tcBU[j].at_i)))
 497         tcr->tcBU[j].bPrint=FALSE;
 498     }
 499
 500   tcr->nQ  = nQ;
 501   tcr->nLJ = nLJ;
 502   tcr->nBU = nBU;
 503   tcr->nIP = nIP;
 504
 505   sfree(inp);
 506
 507   done_warning(wi,FARGS);
 508 }
 509