src/kernel/gctio.c

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