src/gromacs/gmxlib/inputrec.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   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 4.5
  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-2010, 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  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40
  41 #include "typedefs.h"
  42 #include "macros.h"
  43 #include "inputrec.h"
  44 #include "gmx_fatal.h"
  45
  46
  47 /* The minimum number of integration steps required for reasonably accurate
  48  * integration of first and second order coupling algorithms.
  49  */
  50 const int nstmin_berendsen_tcouple =  5;
  51 const int nstmin_berendsen_pcouple = 10;
  52 const int nstmin_harmonic          = 20;
  53
  54 static int nst_wanted(const t_inputrec *ir)
  55 {
  56     if (ir->nstlist > 0)
  57     {
  58         return ir->nstlist;
  59     }
  60     else
  61     {
  62         return 10;
  63     }
  64 }
  65
  66 int ir_optimal_nstcalcenergy(const t_inputrec *ir)
  67 {
  68     return nst_wanted(ir);
  69 }
  70
  71 int tcouple_min_integration_steps(int etc)
  72 {
  73     int n;
  74
  75     switch (etc)
  76     {
  77     case etcNO:
  78         n = 0;
  79         break;
  80     case etcBERENDSEN:
  81     case etcYES:
  82         n = nstmin_berendsen_tcouple;
  83         break;
  84     case etcVRESCALE:
  85         /* V-rescale supports instantaneous rescaling */
  86         n = 0;
  87         break;
  88     case etcNOSEHOOVER:
  89         n = nstmin_harmonic;
  90         break;
  91     case etcANDERSEN:
  92     case etcANDERSENMASSIVE:
  93         n = 1;
  94         break;
  95     default:
  96         gmx_incons("Unknown etc value");
  97         n = 0;
  98     }
  99
 100     return n;
 101 }
 102
 103 int ir_optimal_nsttcouple(const t_inputrec *ir)
 104 {
 105     int  nmin,nwanted,n;
 106     real tau_min;
 107     int  g;
 108
 109     nmin = tcouple_min_integration_steps(ir->etc);
 110
 111     nwanted = nst_wanted(ir);
 112
 113     tau_min = 1e20;
 114     if (ir->etc != etcNO)
 115     {
 116         for(g=0; g<ir->opts.ngtc; g++)
 117         {
 118             if (ir->opts.tau_t[g] > 0)
 119             {
 120                 tau_min = min(tau_min,ir->opts.tau_t[g]);
 121             }
 122         }
 123     }
 124
 125     if (nmin == 0 || ir->delta_t*nwanted <= tau_min)
 126     {
 127         n = nwanted;
 128     }
 129     else
 130     {
 131         n = (int)(tau_min/(ir->delta_t*nmin) + 0.001);
 132         if (n < 1)
 133         {
 134             n = 1;
 135         }
 136         while (nwanted % n != 0)
 137         {
 138             n--;
 139         }
 140     }
 141
 142     return n;
 143 }
 144
 145 int pcouple_min_integration_steps(int epc)
 146 {
 147     int n;
 148
 149     switch (epc)
 150     {
 151     case epcNO:
 152         n = 0;
 153         break;
 154     case etcBERENDSEN:
 155     case epcISOTROPIC:
 156         n = nstmin_berendsen_pcouple;
 157         break;
 158     case epcPARRINELLORAHMAN:
 159     case epcMTTK:
 160         n = nstmin_harmonic;
 161         break;
 162     default:
 163         gmx_incons("Unknown epc value");
 164         n = 0;
 165     }
 166
 167     return n;
 168 }
 169
 170 int ir_optimal_nstpcouple(const t_inputrec *ir)
 171 {
 172     int  nmin,nwanted,n;
 173
 174     nmin = pcouple_min_integration_steps(ir->epc);
 175
 176     nwanted = nst_wanted(ir);
 177
 178     if (nmin == 0 || ir->delta_t*nwanted <= ir->tau_p)
 179     {
 180         n = nwanted;
 181     }
 182     else
 183     {
 184         n = (int)(ir->tau_p/(ir->delta_t*nmin) + 0.001);
 185         if (n < 1)
 186         {
 187             n = 1;
 188         }
 189         while (nwanted % n != 0)
 190         {
 191             n--;
 192         }
 193     }
 194
 195     return n;
 196 }