src/kernel/gctio.c

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