<li><A HREF="#run"><b>run control</b></A> (integrator, tinit, dt, nsteps, init-step, comm-mode, nstcomm, comm-grps)
<li><A HREF="#ld"><b>langevin dynamics</b></A> (bd-fric, ld-seed)
<li><A HREF="#em"><b>energy minimization</b></A> (emtol, emstep, nstcgsteep)
-<li><a HREF="#xmdrun"><b>shell molecular dynamics</b></a>(emtol,niter,fcstep)
+<li><a HREF="#shellmd"><b>shell molecular dynamics</b></a>(emtol,niter,fcstep)
<li><a HREF="#tpi"><b>test particle insertion</b></a>(rtpi)
<li><A HREF="#out"><b>output control</b></A> (nstxout, nstvout, nstfout, nstlog, nstcalcenergy, nstenergy, nstxtcout, xtc-precision, xtc-grps, energygrps)
<li><A HREF="#nl"><b>neighbor searching</b></A> (cutoff-scheme, nstlist, nstcalclr, ns-type, pbc, periodic-molecules, verlet-buffer-drift, rlist, rlistlong)
number is (at least theoretically) more accurate, but slower.</dd>
</dl>
-<A NAME="xmdrun"><br>
+<A NAME="shellmd"><br>
<hr>
<h3>Shell Molecular Dynamics<!--QuietIdx-->shell molecular dynamics<!--EQuietIdx--></h3>
When shells or
<A HREF="#ef">E-yt</A><br>
<A HREF="#ef">E-z</A><br>
<A HREF="#ef">E-zt </A><br>
-<A HREF="#xmdrun">fcstep</A><br>
+<A HREF="#shellmd">fcstep</A><br>
<A HREF="#ewald">fourier-nx</A><br>
<A HREF="#ewald">fourier-ny</A><br>
<A HREF="#ewald">fourier-nz</A><br>
<A HREF="#expanded">mininum-var-min</A><br>
<A HREF="#bond2">morse</A><br>
<A HREF="#em">nbfgscorr</A><br>
-<A HREF="#xmdrun">niter</A><br>
+<A HREF="#shellmd">niter</A><br>
<A HREF="#tc">nh-chain-length</A><br>
<A HREF="#em">nstcgsteep</A><br>
<A HREF="#out">nstcalcenergy</A><br>
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.3.99_development_20071104
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2006, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Groningen Machine for Chemical Simulation
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <string.h>
-#include "typedefs.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "macros.h"
-#include "txtdump.h"
-#include "tpxio.h"
-#include "enxio.h"
-#include "names.h"
-#include "statutil.h"
-#include "copyrite.h"
-#include "random.h"
-
-real ener(matrix P,real e,real e0,int nmol,real kp,real ke,gmx_bool bPScal)
-{
- if (bPScal)
- return (kp*(sqr(P[XX][XX]+P[YY][YY]+P[ZZ][ZZ]-3))+
- ke*sqr(e/nmol-e0));
- else
- return (kp*(sqr(P[XX][XX]-1)+sqr(P[YY][YY]-1)+sqr(P[ZZ][ZZ]-1)+
- sqr(P[XX][YY])+sqr(P[XX][ZZ])+sqr(P[YY][ZZ])) +
- ke*sqr(e/nmol-e0));
-}
-
-void do_sim(char *enx,
- t_topology *top,rvec *x,rvec *v,t_inputrec *ir,matrix box)
-{
- char *tpx = "optwat.tpr";
- char buf[128];
-
- write_tpx(tpx,0,0.0,0.0,ir,box,top->atoms.nr,x,v,NULL,top);
- sprintf(buf,"xmdrun -s %s -e %s >& /dev/null",tpx,enx);
- system(buf);
-}
-
-void get_results(char *enx,real P[],real *epot,int pindex,int eindex)
-{
- ener_file_t fp_ene;
- int nre,step,ndr,i;
- real t;
- gmx_enxnm_t *nms=NULL;
- t_enxframe fr;
-
- fp_ene = open_enx(enx,"r");
-
- do_enxnms(fp_ene,&nre,&nms);
-
- /* Read until the last frame */
- while (do_enx(fp_ene,&fr));
-
- close_enx(fp_ene);
-
- *epot = fr.ener[eindex].e;
- for(i=pindex; (i<pindex+9); i++)
- P[i-pindex] = fr.ener[i].e;
-
- sfree(ener);
-}
-
-void copy_iparams(int nr,t_iparams dest[],t_iparams src[])
-{
- memcpy(dest,src,nr*sizeof(dest[0]));
-}
-
-void rand_step(FILE *fp,int nr,t_iparams ip[],int *seed,real frac)
-{
- int i;
- real ff;
-
- do {
- i = (int) (rando(seed)*nr);
- } while (ip[i].lj.c12 == 0.0);
-
- do {
- ff = frac*rando(seed);
- } while (ff == 0.0);
-
- if (rando(seed) > 0.5) {
- ip[i].lj.c12 *= 1.0+ff;
- fprintf(fp,"Increasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
- }
- else {
- ip[i].lj.c12 *= 1.0-ff;
- fprintf(fp,"Decreasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
- }
-}
-
-void pr_progress(FILE *fp,int nit,tensor P,real epot,real eFF,
- double mc_crit,gmx_bool bConv,gmx_bool bAccept)
-{
- fprintf(fp,"Iter %3d, eFF = %g, Converged = %s, Accepted = %s\n",
- nit,eFF,yesno_names[bConv],yesno_names[bAccept]);
- fprintf(fp,"Epot = %g Pscal = %g, mc_crit = %g\n",epot,
- trace(P)/3,mc_crit);
- pr_rvecs(fp,0,"Pres",P,DIM);
- fprintf(fp,"-----------------------------------------------------\n");
- fflush(fp);
-}
-
-int main(int argc,char *argv[])
-{
- static char *desc[] = {
- "[TT]optwat[tt] optimizes the force field parameter set of a molecular crystal",
- "to reproduce the pressure tensor and experimental energy.[PAR]",
- "Note that for good results the [TT].tpx[tt] file must contain input for a",
- "simulated annealing run, or a single point energy calculation at 0 K"
- };
- t_filenm fnm[] = {
- { efTPX, NULL, NULL, ffREAD },
- { efEDR, "-e", NULL, ffRW },
- { efLOG, "-g", NULL, ffWRITE }
- };
-#define NFILE asize(fnm)
-
- static real epot0 = -57, tol = 1, kT = 0.0;
- static real kp = 1, ke = 100, frac = 0.1;
- static int maxnit = 100, eindex = 5, pindex = 19;
- static int seed = 1993;
- static gmx_bool bPScal = FALSE;
- static t_pargs pa[] = {
- { "-epot0", FALSE, etREAL, {&epot0},
- "Potential energy in kJ/mol" },
- { "-tol", FALSE, etREAL, {&tol},
- "Tolerance for converging" },
- { "-nit", FALSE, etINT, {&maxnit},
- "Max number of iterations" },
- { "-seed", FALSE, etINT, {&seed},
- "Random seed for MC steps" },
- { "-frac", FALSE, etREAL, {&frac},
- "Maximum fraction by which to change parameters. Actual fraction is random between 0 and this parameter" },
- { "-pindex", FALSE, etINT, {&pindex},
- "Index of P[X][X] in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
- { "-eindex", FALSE, etINT, {&pindex},
- "Index of Epot in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
- { "-kp", FALSE, etREAL, {&kp},
- "Force constant for pressure components"},
- { "-ke", FALSE, etREAL, {&ke},
- "Force constant for energy component"},
- { "-kT", FALSE, etREAL, {&kT},
- "Boltzmann Energy for Monte Carlo" },
- { "-pscal", FALSE, etBOOL, {&bPScal},
- "Optimize params for scalar pressure, instead of tensor" }
- };
-
- FILE *fp;
- t_topology top;
- t_tpxheader sh;
- t_inputrec ir;
- t_iparams *ip[2];
- int cur=0;
-#define next (1-cur)
- rvec *xx,*vv;
- matrix box;
- int i,step,natoms,nmol,nit,atnr2;
- real t,lambda,epot,eFF[2];
- double mc_crit=0;
- gmx_bool bConverged,bAccept;
- tensor P;
-
- CopyRight(stdout,argv[0]);
- parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
- asize(desc),desc,0,NULL);
-
- /* Read initial topology and coordaintes etc. */
- read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&sh,TRUE,NULL,NULL);
- snew(xx,sh.natoms);
- snew(vv,sh.natoms);
- read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,box,&natoms,
- xx,vv,NULL,&top);
-
- /* Open log file and print options */
- fp = ftp2FILE(efLOG,NFILE,fnm,"w");
- fprintf(fp,"%s started with the following parameters\n",argv[0]);
- fprintf(fp,"epot = %8g ke = %8g kp = %8g\n",epot0,ke,kp);
- fprintf(fp,"maxnit = %8d tol = %8g seed = %8d\n",maxnit,tol,seed);
- fprintf(fp,"frac = %8g pindex = %8d eindex = %8d\n",frac,pindex,eindex);
- fprintf(fp,"kT = %8g pscal = %8s\n",kT,gmx_bool_names[bPScal]);
-
- /* Unpack some topology numbers */
- nmol = top.blocks[ebMOLS].nr;
- atnr2 = top.idef.atnr*top.idef.atnr;
-
- /* Copy input params */
- snew(ip[cur],atnr2);
- snew(ip[next],atnr2);
- copy_iparams(atnr2,ip[cur],top.idef.iparams);
- copy_iparams(atnr2,ip[next],top.idef.iparams);
-
- /* Loop over iterations */
- nit = 0;
- do {
- if (nit > 0) {
- /* Do random step */
- rand_step(fp,atnr2,ip[next],&seed,frac);
- copy_iparams(atnr2,top.idef.iparams,ip[next]);
- }
- do_sim(ftp2fn(efEDR,NFILE,fnm),&top,xx,vv,&ir,box);
-
- get_results(ftp2fn(efEDR,NFILE,fnm),P[0],&epot,pindex,eindex);
-
- /* Calculate penalty */
- eFF[(nit > 0) ? next : cur] = ener(P,epot,epot0,nmol,kp,ke,bPScal);
-
- bConverged = (eFF[(nit > 0) ? next : cur] < tol);
-
- if (nit > 0) {
- /* Do Metropolis criterium */
- if (kT > 0)
- mc_crit = exp(-(eFF[next]-eFF[cur])/kT);
- bAccept = ((eFF[next] < eFF[cur]) ||
- ((kT > 0) && (mc_crit > rando(&seed))));
- pr_progress(fp,nit,P,epot/nmol,eFF[next],mc_crit,
- bConverged,bAccept);
- if (bAccept) {
- /* Better params! */
- cur = next;
- }
- else {
- /* Restore old parameters */
- copy_iparams(atnr2,ip[next],ip[cur]);
- }
- }
- else
- pr_progress(fp,nit,P,epot/nmol,eFF[cur],mc_crit,bConverged,FALSE);
-
- nit++;
- } while (!bConverged && (nit < maxnit));
-
- for(i=0; (i<atnr2); i++)
- pr_iparams(fp,F_LJ,&ip[cur][i]);
-
- ffclose(fp);
-
- gmx_thanx(stderr);
-
- return 0;
-}
{ efXVG, "-tpid", "tpidist", ffOPTWR },
{ efEDI, "-ei", "sam", ffOPTRD },
{ efXVG, "-eo", "edsam", ffOPTWR },
- { efGCT, "-j", "wham", ffOPTRD },
- { efGCT, "-jo", "bam", ffOPTWR },
- { efXVG, "-ffout", "gct", ffOPTWR },
{ efXVG, "-devout", "deviatie", ffOPTWR },
{ efXVG, "-runav", "runaver", ffOPTWR },
{ efXVG, "-px", "pullx", ffOPTWR },
{ efXVG, "-bfield", "benchfld", ffOPTWR },
{ efXVG, "-btpi", "benchtpi", ffOPTWR },
{ efXVG, "-btpid", "benchtpid", ffOPTWR },
- { efGCT, "-bjo", "bench", ffOPTWR },
- { efXVG, "-bffout", "benchgct", ffOPTWR },
{ efXVG, "-bdevout", "benchdev", ffOPTWR },
{ efXVG, "-brunav", "benchrnav", ffOPTWR },
{ efXVG, "-bpx", "benchpx", ffOPTWR },
deffile[efNR] =
{
{ eftASC, ".mdp", "grompp", "-f", "grompp input file with MD parameters" },
- { eftASC, ".gct", "gct", "-f", "General coupling stuff"},
{ eftGEN, ".???", "traj", "-f",
"Trajectory: xtc trr trj gro g96 pdb cpt", NTRXS, trxs },
{ eftGEN, ".???", "trajout", "-f",
/* this enum should correspond to the array deffile in gmxlib/filenm.c */
enum {
- efMDP, efGCT,
+ efMDP,
efTRX, efTRO, efTRN, efTRR, efTRJ, efXTC, efG87,
efEDR,
efSTX, efSTO, efGRO, efG96, efPDB, efBRK, efENT, efESP, efPQR, efXYZ,
/* Energy group pair flags */
int *egp_flags;
- /* xmdrun flexible constraints */
+ /* Shell molecular dynamics flexible constraints */
real fc_stepsize;
/* Generalized born implicit solvent */
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "typedefs.h"
-#include "xmdrun.h"
-#include "futil.h"
-#include "xvgr.h"
-#include "macros.h"
-#include "physics.h"
-#include "network.h"
-#include "smalloc.h"
-#include "mdrun.h"
-#include "string2.h"
-#include "readinp.h"
-#include "filenm.h"
-#include "update.h"
-#include "vec.h"
-#include "main.h"
-#include "txtdump.h"
-
-/*#define DEBUGGCT*/
-t_coupl_rec *init_coupling(FILE *log, int nfile, const t_filenm fnm[],
- t_commrec *cr, t_forcerec *fr,
- t_mdatoms *md, t_idef *idef)
-{
- int i, nc, index, j;
- int ati, atj;
- t_coupl_rec *tcr;
-
- snew(tcr, 1);
- if (MASTER(cr))
- {
- read_gct (opt2fn("-j", nfile, fnm), tcr);
- write_gct(opt2fn("-jo", nfile, fnm), tcr, idef);
- }
- /* Update all processors with coupling info */
- if (PAR(cr))
- {
- comm_tcr(log, cr, &tcr);
- }
-
- /* Copy information from the coupling to the force field stuff */
- copy_ff(tcr, fr, md, idef);
-
- return tcr;
-}
-
-static real Ecouple(t_coupl_rec *tcr, real ener[])
-{
- if (tcr->bInter)
- {
- return ener[F_COUL_SR]+ener[F_LJ]+ener[F_COUL_LR]+ener[F_LJ_LR]+
- ener[F_RF_EXCL]+ener[F_COUL_RECIP];
- }
- else
- {
- return ener[F_EPOT];
- }
-}
-
-static char *mk_gct_nm(const char *fn, int ftp, int ati, int atj)
-{
- static char buf[256];
-
- strcpy(buf, fn);
- if (atj == -1)
- {
- sprintf(buf+strlen(fn)-4, "%d.%s", ati, ftp2ext(ftp));
- }
- else
- {
- sprintf(buf+strlen(fn)-4, "%d_%d.%s", ati, atj, ftp2ext(ftp));
- }
-
- return buf;
-}
-
-static void pr_ff(t_coupl_rec *tcr, real time, t_idef *idef,
- t_commrec *cr, int nfile, const t_filenm fnm[],
- const output_env_t oenv)
-{
- static FILE *prop;
- static FILE **out = NULL;
- static FILE **qq = NULL;
- static FILE **ip = NULL;
- t_coupl_LJ *tclj;
- t_coupl_BU *tcbu;
- char buf[256];
- const char *leg[] = { "C12", "C6" };
- const char *eleg[] = { "Epsilon", "Sigma" };
- const char *bleg[] = { "A", "B", "C" };
- char **raleg;
- int i, j, index;
-
- if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL))
- {
- prop = xvgropen(opt2fn("-runav", nfile, fnm),
- "Properties and Running Averages", "Time (ps)", "", oenv);
- snew(raleg, 2*eoObsNR);
- for (i = j = 0; (i < eoObsNR); i++)
- {
- if (tcr->bObsUsed[i])
- {
- raleg[j++] = strdup(eoNames[i]);
- sprintf(buf, "RA-%s", eoNames[i]);
- raleg[j++] = strdup(buf);
- }
- }
- xvgr_legend(prop, j, (const char**)raleg, oenv);
- for (i = 0; (i < j); i++)
- {
- sfree(raleg[i]);
- }
- sfree(raleg);
-
- if (tcr->nLJ)
- {
- snew(out, tcr->nLJ);
- for (i = 0; (i < tcr->nLJ); i++)
- {
- if (tcr->tcLJ[i].bPrint)
- {
- tclj = &(tcr->tcLJ[i]);
- out[i] =
- xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm),
- efXVG, tclj->at_i, tclj->at_j),
- "General Coupling Lennard Jones", "Time (ps)",
- "Force constant (units)", oenv);
- fprintf(out[i], "@ subtitle \"Interaction between types %d and %d\"\n",
- tclj->at_i, tclj->at_j);
- if (tcr->combrule == 1)
- {
- xvgr_legend(out[i], asize(leg), leg, oenv);
- }
- else
- {
- xvgr_legend(out[i], asize(eleg), eleg, oenv);
- }
- fflush(out[i]);
- }
- }
- }
- else if (tcr->nBU)
- {
- snew(out, tcr->nBU);
- for (i = 0; (i < tcr->nBU); i++)
- {
- if (tcr->tcBU[i].bPrint)
- {
- tcbu = &(tcr->tcBU[i]);
- out[i] =
- xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG,
- tcbu->at_i, tcbu->at_j),
- "General Coupling Buckingham", "Time (ps)",
- "Force constant (units)", oenv);
- fprintf(out[i], "@ subtitle \"Interaction between types %d and %d\"\n",
- tcbu->at_i, tcbu->at_j);
- xvgr_legend(out[i], asize(bleg), bleg, oenv);
- fflush(out[i]);
- }
- }
- }
- snew(qq, tcr->nQ);
- for (i = 0; (i < tcr->nQ); i++)
- {
- if (tcr->tcQ[i].bPrint)
- {
- qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG,
- tcr->tcQ[i].at_i, -1),
- "General Coupling Charge", "Time (ps)", "Charge (e)",
- oenv);
- fprintf(qq[i], "@ subtitle \"Type %d\"\n", tcr->tcQ[i].at_i);
- fflush(qq[i]);
- }
- }
- snew(ip, tcr->nIP);
- for (i = 0; (i < tcr->nIP); i++)
- {
- sprintf(buf, "gctIP%d", tcr->tIP[i].type);
- ip[i] = xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG, 0, -1),
- "General Coupling iparams", "Time (ps)", "ip ()", oenv);
- index = tcr->tIP[i].type;
- fprintf(ip[i], "@ subtitle \"Coupling to %s\"\n",
- interaction_function[idef->functype[index]].longname);
- fflush(ip[i]);
- }
- }
- /* Write properties to file */
- fprintf(prop, "%10.3f", time);
- for (i = 0; (i < eoObsNR); i++)
- {
- if (tcr->bObsUsed[i])
- {
- fprintf(prop, " %10.3e %10.3e", tcr->act_value[i], tcr->av_value[i]);
- }
- }
- fprintf(prop, "\n");
- fflush(prop);
-
- for (i = 0; (i < tcr->nLJ); i++)
- {
- tclj = &(tcr->tcLJ[i]);
- if (tclj->bPrint)
- {
- if (tcr->combrule == 1)
- {
- fprintf(out[i], "%14.7e %14.7e %14.7e\n",
- time, tclj->c12, tclj->c6);
- }
- else
- {
- double sigma = pow(tclj->c12/tclj->c6, 1.0/6.0);
- double epsilon = 0.25*sqr(tclj->c6)/tclj->c12;
- fprintf(out[i], "%14.7e %14.7e %14.7e\n",
- time, epsilon, sigma);
- }
- fflush(out[i]);
- }
- }
- for (i = 0; (i < tcr->nBU); i++)
- {
- tcbu = &(tcr->tcBU[i]);
- if (tcbu->bPrint)
- {
- fprintf(out[i], "%14.7e %14.7e %14.7e %14.7e\n",
- time, tcbu->a, tcbu->b, tcbu->c);
- fflush(out[i]);
- }
- }
- for (i = 0; (i < tcr->nQ); i++)
- {
- if (tcr->tcQ[i].bPrint)
- {
- fprintf(qq[i], "%14.7e %14.7e\n", time, tcr->tcQ[i].Q);
- fflush(qq[i]);
- }
- }
- for (i = 0; (i < tcr->nIP); i++)
- {
- fprintf(ip[i], "%10g ", time);
- index = tcr->tIP[i].type;
- switch (idef->functype[index])
- {
- case F_BONDS:
- fprintf(ip[i], "%10g %10g\n", tcr->tIP[i].iprint.harmonic.krA,
- tcr->tIP[i].iprint.harmonic.rA);
- break;
- default:
- break;
- }
- fflush(ip[i]);
- }
-}
-
-static void pr_dev(t_coupl_rec *tcr,
- real t, real dev[eoObsNR], t_commrec *cr, int nfile,
- const t_filenm fnm[], const output_env_t oenv)
-{
- static FILE *fp = NULL;
- char **ptr;
- int i, j;
-
- if (!fp)
- {
- fp = xvgropen(opt2fn("-devout", nfile, fnm),
- "Deviations from target value", "Time (ps)", "", oenv);
- snew(ptr, eoObsNR);
- for (i = j = 0; (i < eoObsNR); i++)
- {
- if (tcr->bObsUsed[i])
- {
- ptr[j++] = strdup(eoNames[i]);
- }
- }
- xvgr_legend(fp, j, (const char**)ptr, oenv);
- for (i = 0; i < j; i++)
- {
- sfree(ptr[i]);
- }
- sfree(ptr);
- }
- fprintf(fp, "%10.3f", t);
- for (i = 0; (i < eoObsNR); i++)
- {
- if (tcr->bObsUsed[i])
- {
- fprintf(fp, " %10.3e", dev[i]);
- }
- }
- fprintf(fp, "\n");
- fflush(fp);
-}
-
-static void upd_nbfplj(FILE *log, real *nbfp, int atnr, real f6[], real f12[],
- int combrule)
-{
- double *sigma, *epsilon, c6, c12, eps, sig, sig6;
- int n, m, k;
-
- /* Update the nonbonded force parameters */
- switch (combrule)
- {
- case 1:
- for (k = n = 0; (n < atnr); n++)
- {
- for (m = 0; (m < atnr); m++, k++)
- {
- C6 (nbfp, atnr, n, m) *= f6[k];
- C12(nbfp, atnr, n, m) *= f12[k];
- }
- }
- break;
- case 2:
- case 3:
- /* Convert to sigma and epsilon */
- snew(sigma, atnr);
- snew(epsilon, atnr);
- for (n = 0; (n < atnr); n++)
- {
- k = n*(atnr+1);
- c6 = C6 (nbfp, atnr, n, n) * f6[k];
- c12 = C12(nbfp, atnr, n, n) * f12[k];
- if ((c6 == 0) || (c12 == 0))
- {
- gmx_fatal(FARGS, "You can not use combination rule %d with zero C6 (%f) or C12 (%f)", combrule, c6, c12);
- }
- sigma[n] = pow(c12/c6, 1.0/6.0);
- epsilon[n] = 0.25*(c6*c6/c12);
- }
- for (k = n = 0; (n < atnr); n++)
- {
- for (m = 0; (m < atnr); m++, k++)
- {
- eps = sqrt(epsilon[n]*epsilon[m]);
- if (combrule == 2)
- {
- sig = 0.5*(sigma[n]+sigma[m]);
- }
- else
- {
- sig = sqrt(sigma[n]*sigma[m]);
- }
- sig6 = pow(sig, 6.0);
- /* nbfp now includes the 6.0/12.0 derivative prefactors */
- C6 (nbfp, atnr, n, m) = 4*eps*sig6/6.0;
- C12(nbfp, atnr, n, m) = 4*eps*sig6*sig6/12.0;
- }
- }
- sfree(sigma);
- sfree(epsilon);
- break;
- default:
- gmx_fatal(FARGS, "Combination rule should be 1,2 or 3 instead of %d",
- combrule);
- }
-}
-
-static void upd_nbfpbu(FILE *log, real *nbfp, int atnr,
- real fa[], real fb[], real fc[])
-{
- int n, m, k;
-
- /* Update the nonbonded force parameters */
- for (k = n = 0; (n < atnr); n++)
- {
- for (m = 0; (m < atnr); m++, k++)
- {
- BHAMA(nbfp, atnr, n, m) *= fa[k];
- BHAMB(nbfp, atnr, n, m) *= fb[k];
- BHAMC(nbfp, atnr, n, m) *= fc[k];
- }
- }
-}
-
-void gprod(t_commrec *cr, int n, real f[])
-{
- /* Compute the global product of all elements in an array
- * such that after gprod f[i] = PROD_j=1,nprocs f[i][j]
- */
- static int nbuf = 0;
- static real *buf = NULL;
- int i;
-
- if (n > nbuf)
- {
- nbuf = n;
- srenew(buf, nbuf);
- }
-
-#ifdef GMX_MPI
-#ifdef GMX_DOUBLE
- MPI_Allreduce(f, buf, n, MPI_DOUBLE, MPI_PROD, cr->mpi_comm_mygroup);
-#else
- MPI_Allreduce(f, buf, n, MPI_FLOAT, MPI_PROD, cr->mpi_comm_mygroup);
-#endif
- for (i = 0; (i < n); i++)
- {
- f[i] = buf[i];
- }
-#endif
-}
-
-static void set_factor_matrix(int ntypes, real f[], real fmult, int ati, int atj)
-{
-#define FMIN 0.95
-#define FMAX 1.05
- int i;
-
- fmult = min(FMAX, max(FMIN, fmult));
- if (atj != -1)
- {
- f[ntypes*ati+atj] *= fmult;
- f[ntypes*atj+ati] *= fmult;
- }
- else
- {
- for (i = 0; (i < ntypes); i++)
- {
- f[ntypes*ati+i] *= fmult;
- f[ntypes*i+ati] *= fmult;
- }
- }
-#undef FMIN
-#undef FMAX
-}
-
-static real calc_deviation(real xav, real xt, real x0)
-{
- /* This may prevent overshooting in GCT coupling... */
-
- /* real dev;
-
- if (xav > x0) {
- if (xt > x0)
- dev = min(xav-x0,xt-x0);
- else
- dev = 0;
- }
- else {
- if (xt < x0)
- dev = max(xav-x0,xt-x0);
- else
- dev = 0;
- }
- */
- return x0-xav;
-}
-
-static real calc_dist(FILE *log, rvec x[])
-{
- static gmx_bool bFirst = TRUE;
- static gmx_bool bDist;
- static int i1, i2;
- char *buf;
- rvec dx;
-
- if (bFirst)
- {
- if ((buf = getenv("DISTGCT")) == NULL)
- {
- bDist = FALSE;
- }
- else
- {
- bDist = (sscanf(buf, "%d%d", &i1, &i2) == 2);
- if (bDist)
- {
- fprintf(log, "GCT: Will couple to distance between %d and %d\n", i1, i2);
- }
- else
- {
- fprintf(log, "GCT: Will not couple to distances\n");
- }
- }
- bFirst = FALSE;
- }
- if (bDist)
- {
- rvec_sub(x[i1], x[i2], dx);
- return norm(dx);
- }
- else
- {
- return 0.0;
- }
-}
-
-real run_aver(real old, real cur, int step, int nmem)
-{
- nmem = max(1, nmem);
-
- return ((nmem-1)*old+cur)/nmem;
-}
-
-static void set_act_value(t_coupl_rec *tcr, int index, real val, int step)
-{
- tcr->act_value[index] = val;
- tcr->av_value[index] = run_aver(tcr->av_value[index], val, step, tcr->nmemory);
-}
-
-static void upd_f_value(FILE *log, int atnr, real xi, real dt, real factor,
- real ff[], int ati, int atj)
-{
- real fff;
-
- if (xi != 0)
- {
- fff = 1 + (dt/xi) * factor;
- if (fff > 0)
- {
- set_factor_matrix(atnr, ff, sqrt(fff), ati, atj);
- }
- }
-}
-
-static void dump_fm(FILE *fp, int n, real f[], char *s)
-{
- int i, j;
-
- fprintf(fp, "Factor matrix (all numbers -1) %s\n", s);
- for (i = 0; (i < n); i++)
- {
- for (j = 0; (j < n); j++)
- {
- fprintf(fp, " %10.3e", f[n*i+j]-1.0);
- }
- fprintf(fp, "\n");
- }
-}
-
-void do_coupling(FILE *log, const output_env_t oenv, int nfile,
- const t_filenm fnm[], t_coupl_rec *tcr, real t,
- int step, real ener[], t_forcerec *fr, t_inputrec *ir,
- gmx_bool bMaster, t_mdatoms *md, t_idef *idef, real mu_aver, int nmols,
- t_commrec *cr, matrix box, tensor virial,
- tensor pres, rvec mu_tot,
- rvec x[], rvec f[], gmx_bool bDoIt)
-{
-#define enm2Debye 48.0321
-#define d2e(x) (x)/enm2Debye
-#define enm2kjmol(x) (x)*0.0143952 /* = 2.0*4.0*M_PI*EPSILON0 */
-
- static real *f6, *f12, *fa, *fb, *fc, *fq;
- static gmx_bool bFirst = TRUE;
-
- int i, j, ati, atj, atnr2, type, ftype;
- real deviation[eoObsNR], prdev[eoObsNR], epot0, dist, rmsf;
- real ff6, ff12, ffa, ffb, ffc, ffq, factor, dt, mu_ind;
- real Epol, Eintern, Virial, muabs, xiH = -1, xiS = -1, xi6, xi12;
- rvec fmol[2];
- gmx_bool bTest, bPrint;
- t_coupl_LJ *tclj;
- t_coupl_BU *tcbu;
- t_coupl_Q *tcq;
- t_coupl_iparams *tip;
-
- atnr2 = idef->atnr * idef->atnr;
- if (bFirst)
- {
- if (PAR(cr))
- {
- fprintf(log, "GCT: this is parallel\n");
- }
- else
- {
- fprintf(log, "GCT: this is not parallel\n");
- }
- fflush(log);
- snew(f6, atnr2);
- snew(f12, atnr2);
- snew(fa, atnr2);
- snew(fb, atnr2);
- snew(fc, atnr2);
- snew(fq, idef->atnr);
-
- if (tcr->bVirial)
- {
- int nrdf = 0;
- real TTT = 0;
- real Vol = det(box);
-
- for (i = 0; (i < ir->opts.ngtc); i++)
- {
- nrdf += ir->opts.nrdf[i];
- TTT += ir->opts.nrdf[i]*ir->opts.ref_t[i];
- }
- TTT /= nrdf;
-
- /* Calculate reference virial from reference temperature and pressure */
- tcr->ref_value[eoVir] = 0.5*BOLTZ*nrdf*TTT - (3.0/2.0)*
- Vol*tcr->ref_value[eoPres];
-
- fprintf(log, "GCT: TTT = %g, nrdf = %d, vir0 = %g, Vol = %g\n",
- TTT, nrdf, tcr->ref_value[eoVir], Vol);
- fflush(log);
- }
- bFirst = FALSE;
- }
-
- bPrint = MASTER(cr) && do_per_step(step, ir->nstlog);
- dt = ir->delta_t;
-
- /* Initiate coupling to the reference pressure and temperature to start
- * coupling slowly.
- */
- if (step == 0)
- {
- for (i = 0; (i < eoObsNR); i++)
- {
- tcr->av_value[i] = tcr->ref_value[i];
- }
- if ((tcr->ref_value[eoDipole]) != 0.0)
- {
- mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
- Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
- tcr->av_value[eoEpot] -= Epol;
- fprintf(log, "GCT: mu_aver = %g(D), mu_ind = %g(D), Epol = %g (kJ/mol)\n",
- mu_aver*enm2Debye, mu_ind*enm2Debye, Epol);
- }
- }
-
- /* We want to optimize the LJ params, usually to the Vaporization energy
- * therefore we only count intermolecular degrees of freedom.
- * Note that this is now optional. switch UseEinter to yes in your gct file
- * if you want this.
- */
- dist = calc_dist(log, x);
- muabs = norm(mu_tot);
- Eintern = Ecouple(tcr, ener)/nmols;
- Virial = virial[XX][XX]+virial[YY][YY]+virial[ZZ][ZZ];
-
- /*calc_force(md->nr,f,fmol);*/
- clear_rvec(fmol[0]);
-
- /* Use a memory of tcr->nmemory steps, so we actually couple to the
- * average observable over the last tcr->nmemory steps. This may help
- * in avoiding local minima in parameter space.
- */
- set_act_value(tcr, eoPres, ener[F_PRES], step);
- set_act_value(tcr, eoEpot, Eintern, step);
- set_act_value(tcr, eoVir, Virial, step);
- set_act_value(tcr, eoDist, dist, step);
- set_act_value(tcr, eoMu, muabs, step);
- set_act_value(tcr, eoFx, fmol[0][XX], step);
- set_act_value(tcr, eoFy, fmol[0][YY], step);
- set_act_value(tcr, eoFz, fmol[0][ZZ], step);
- set_act_value(tcr, eoPx, pres[XX][XX], step);
- set_act_value(tcr, eoPy, pres[YY][YY], step);
- set_act_value(tcr, eoPz, pres[ZZ][ZZ], step);
-
- epot0 = tcr->ref_value[eoEpot];
- /* If dipole != 0.0 assume we want to use polarization corrected coupling */
- if ((tcr->ref_value[eoDipole]) != 0.0)
- {
- mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
-
- Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
-
- epot0 -= Epol;
- if (debug)
- {
- fprintf(debug, "mu_ind: %g (%g D) mu_aver: %g (%g D)\n",
- mu_ind, mu_ind*enm2Debye, mu_aver, mu_aver*enm2Debye);
- fprintf(debug, "Eref %g Epol %g Erunav %g Eact %g\n",
- tcr->ref_value[eoEpot], Epol, tcr->av_value[eoEpot],
- tcr->act_value[eoEpot]);
- }
- }
-
- if (bPrint)
- {
- pr_ff(tcr, t, idef, cr, nfile, fnm, oenv);
- }
- /* Calculate the deviation of average value from the target value */
- for (i = 0; (i < eoObsNR); i++)
- {
- deviation[i] = calc_deviation(tcr->av_value[i], tcr->act_value[i],
- tcr->ref_value[i]);
- prdev[i] = tcr->ref_value[i] - tcr->act_value[i];
- }
- deviation[eoEpot] = calc_deviation(tcr->av_value[eoEpot], tcr->act_value[eoEpot],
- epot0);
- prdev[eoEpot] = epot0 - tcr->act_value[eoEpot];
-
- if (bPrint)
- {
- pr_dev(tcr, t, prdev, cr, nfile, fnm, oenv);
- }
-
- /* First set all factors to 1 */
- for (i = 0; (i < atnr2); i++)
- {
- f6[i] = f12[i] = fa[i] = fb[i] = fc[i] = 1.0;
- }
- for (i = 0; (i < idef->atnr); i++)
- {
- fq[i] = 1.0;
- }
-
- /* Now compute the actual coupling compononents */
- if (!fr->bBHAM)
- {
- if (bDoIt)
- {
- for (i = 0; (i < tcr->nLJ); i++)
- {
- tclj = &(tcr->tcLJ[i]);
-
- ati = tclj->at_i;
- atj = tclj->at_j;
-
- ff6 = ff12 = 1.0;
-
- if (tclj->eObs == eoForce)
- {
- gmx_fatal(FARGS, "Hack code for this to work again ");
- if (debug)
- {
- fprintf(debug, "Have computed derivatives: xiH = %g, xiS = %g\n", xiH, xiS);
- }
- if (ati == 1)
- {
- /* Hydrogen */
- ff12 += xiH;
- }
- else if (ati == 2)
- {
- /* Shell */
- ff12 += xiS;
- }
- else
- {
- gmx_fatal(FARGS, "No H, no Shell, edit code at %s, line %d\n",
- __FILE__, __LINE__);
- }
- if (ff6 > 0)
- {
- set_factor_matrix(idef->atnr, f6, sqrt(ff6), ati, atj);
- }
- if (ff12 > 0)
- {
- set_factor_matrix(idef->atnr, f12, sqrt(ff12), ati, atj);
- }
- }
- else
- {
- if (debug)
- {
- fprintf(debug, "tcr->tcLJ[%d].xi_6 = %g, xi_12 = %g deviation = %g\n", i,
- tclj->xi_6, tclj->xi_12, deviation[tclj->eObs]);
- }
- factor = deviation[tclj->eObs];
-
- upd_f_value(log, idef->atnr, tclj->xi_6, dt, factor, f6, ati, atj);
- upd_f_value(log, idef->atnr, tclj->xi_12, dt, factor, f12, ati, atj);
- }
- }
- }
- if (PAR(cr))
- {
- gprod(cr, atnr2, f6);
- gprod(cr, atnr2, f12);
-#ifdef DEBUGGCT
- dump_fm(log, idef->atnr, f6, "f6");
- dump_fm(log, idef->atnr, f12, "f12");
-#endif
- }
- upd_nbfplj(log, fr->nbfp, idef->atnr, f6, f12, tcr->combrule);
-
- /* Copy for printing */
- for (i = 0; (i < tcr->nLJ); i++)
- {
- tclj = &(tcr->tcLJ[i]);
- ati = tclj->at_i;
- atj = tclj->at_j;
- if (atj == -1)
- {
- atj = ati;
- }
- /* nbfp now includes the 6.0/12.0 derivative prefactors */
- tclj->c6 = C6(fr->nbfp, fr->ntype, ati, atj)/6.0;
- tclj->c12 = C12(fr->nbfp, fr->ntype, ati, atj)/12.0;
- }
- }
- else
- {
- if (bDoIt)
- {
- for (i = 0; (i < tcr->nBU); i++)
- {
- tcbu = &(tcr->tcBU[i]);
- factor = deviation[tcbu->eObs];
- ati = tcbu->at_i;
- atj = tcbu->at_j;
-
- upd_f_value(log, idef->atnr, tcbu->xi_a, dt, factor, fa, ati, atj);
- upd_f_value(log, idef->atnr, tcbu->xi_b, dt, factor, fb, ati, atj);
- upd_f_value(log, idef->atnr, tcbu->xi_c, dt, factor, fc, ati, atj);
- }
- }
- if (PAR(cr))
- {
- gprod(cr, atnr2, fa);
- gprod(cr, atnr2, fb);
- gprod(cr, atnr2, fc);
- }
- upd_nbfpbu(log, fr->nbfp, idef->atnr, fa, fb, fc);
- /* Copy for printing */
- for (i = 0; (i < tcr->nBU); i++)
- {
- tcbu = &(tcr->tcBU[i]);
- ati = tcbu->at_i;
- atj = tcbu->at_j;
- if (atj == -1)
- {
- atj = ati;
- }
- /* nbfp now includes the 6.0 derivative prefactors */
- tcbu->a = BHAMA(fr->nbfp, fr->ntype, ati, atj);
- tcbu->b = BHAMB(fr->nbfp, fr->ntype, ati, atj);
- tcbu->c = BHAMC(fr->nbfp, fr->ntype, ati, atj)/6.0;
- if (debug)
- {
- fprintf(debug, "buck (type=%d) = %e, %e, %e\n",
- tcbu->at_i, tcbu->a, tcbu->b, tcbu->c);
- }
- }
- }
- if (bDoIt)
- {
- for (i = 0; (i < tcr->nQ); i++)
- {
- tcq = &(tcr->tcQ[i]);
- if (tcq->xi_Q)
- {
- ffq = 1.0 + (dt/tcq->xi_Q) * deviation[tcq->eObs];
- }
- else
- {
- ffq = 1.0;
- }
- fq[tcq->at_i] *= ffq;
- }
- }
- if (PAR(cr))
- {
- gprod(cr, idef->atnr, fq);
- }
-
- for (j = 0; (j < md->nr); j++)
- {
- md->chargeA[j] *= fq[md->typeA[j]];
- }
- for (i = 0; (i < tcr->nQ); i++)
- {
- tcq = &(tcr->tcQ[i]);
- for (j = 0; (j < md->nr); j++)
- {
- if (md->typeA[j] == tcq->at_i)
- {
- tcq->Q = md->chargeA[j];
- break;
- }
- }
- if (j == md->nr)
- {
- gmx_fatal(FARGS, "Coupling type %d not found", tcq->at_i);
- }
- }
- for (i = 0; (i < tcr->nIP); i++)
- {
- tip = &(tcr->tIP[i]);
- type = tip->type;
- ftype = idef->functype[type];
- factor = dt*deviation[tip->eObs];
-
- switch (ftype)
- {
- case F_BONDS:
- if (tip->xi.harmonic.krA)
- {
- idef->iparams[type].harmonic.krA *= (1+factor/tip->xi.harmonic.krA);
- }
- if (tip->xi.harmonic.rA)
- {
- idef->iparams[type].harmonic.rA *= (1+factor/tip->xi.harmonic.rA);
- }
- break;
- default:
- break;
- }
- tip->iprint = idef->iparams[type];
- }
-}
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "typedefs.h"
-#include "xmdrun.h"
-#include "futil.h"
-#include "xvgr.h"
-#include "macros.h"
-#include "physics.h"
-#include "network.h"
-#include "smalloc.h"
-#include "string2.h"
-#include "readinp.h"
-#include "filenm.h"
-#include "names.h"
-#include "gmxfio.h"
-
-const char *eoNames[eoNR] = {
- "Pres", "Epot", "Vir", "Dist", "Mu", "Force", "Fx", "Fy", "Fz",
- "Px", "Py", "Pz",
- "Polarizability", "Dipole", "Memory", "UseEinter", "UseVirial",
- "CombinationRule"
-};
-
-static int Name2eo(char *s)
-{
- int i, res;
-
- res = -1;
-
- for (i = 0; (i < eoNR); i++)
- {
- if (gmx_strcasecmp(s, eoNames[i]) == 0)
- {
- res = i;
- fprintf(stderr, "Coupling to observable %d (%s)\n", res, eoNames[res]);
- break;
- }
- }
-
- return res;
-}
-
-#define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
-#define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
-#define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
-
-static void low_comm_tcr(t_commrec *cr, t_coupl_rec *tcr)
-{
- nblock_bc(cr, eoObsNR, tcr->ref_value);
-
- block_bc(cr, tcr->nLJ);
- snew_bc(cr, tcr->tcLJ, tcr->nLJ);
- nblock_bc(cr, tcr->nLJ, tcr->tcLJ);
-
- block_bc(cr, tcr->nBU);
- snew_bc(cr, tcr->tcBU, tcr->nBU);
- nblock_bc(cr, tcr->nBU, tcr->tcBU);
-
- block_bc(cr, tcr->nQ);
- snew_bc(cr, tcr->tcQ, tcr->nQ);
- nblock_bc(cr, tcr->nQ, tcr->tcQ);
-
- block_bc(cr, tcr->nmemory);
- block_bc(cr, tcr->bInter);
- block_bc(cr, tcr->bVirial);
- block_bc(cr, tcr->combrule);
-}
-
-void comm_tcr(FILE *log, t_commrec *cr, t_coupl_rec **tcr)
-{
- if (!MASTER(cr))
- {
- snew(*tcr, 1);
- }
-
- low_comm_tcr(cr, *tcr);
-}
-
-static void clear_lj(t_coupl_LJ *tc)
-{
- tc->at_i = 0;
- tc->at_j = 0;
- tc->eObs = -1;
- tc->bPrint = TRUE;
- tc->c6 = 0.0;
- tc->c12 = 0.0;
- tc->xi_6 = 0.0;
- tc->xi_12 = 0.0;
-}
-
-static void clear_bu(t_coupl_BU *tc)
-{
- tc->at_i = 0;
- tc->at_j = 0;
- tc->eObs = -1;
- tc->bPrint = TRUE;
- tc->a = 0.0;
- tc->b = 0.0;
- tc->c = 0.0;
- tc->xi_a = 0.0;
- tc->xi_b = 0.0;
- tc->xi_c = 0.0;
-}
-
-static void clear_q(t_coupl_Q *tc)
-{
- tc->at_i = 0;
- tc->eObs = -1;
- tc->bPrint = TRUE;
- tc->Q = 0.0;
- tc->xi_Q = 0.0;
-}
-
-void copy_ff(t_coupl_rec *tcr, t_forcerec *fr, t_mdatoms *md, t_idef *idef)
-{
- int i, j, ati, atj, type;
- t_coupl_LJ *tclj;
- t_coupl_BU *tcbu;
- t_coupl_Q *tcq;
-
- /* Save values for printing */
- for (i = 0; (i < tcr->nLJ); i++)
- {
- tclj = &(tcr->tcLJ[i]);
-
- ati = tclj->at_i;
- atj = tclj->at_j;
- if (atj == -1)
- {
- atj = ati;
- }
- /* nbfp now includes the 6.0/12.0 derivative prefactors */
- tclj->c6 = C6(fr->nbfp, fr->ntype, ati, atj)/6.0;
- tclj->c12 = C12(fr->nbfp, fr->ntype, ati, atj)/12.0;
- }
-
- for (i = 0; (i < tcr->nBU); i++)
- {
- tcbu = &(tcr->tcBU[i]);
-
- ati = tcbu->at_i;
- atj = tcbu->at_j;
- if (atj == -1)
- {
- atj = ati;
- }
- /* nbfp now includes the 6.0 derivative prefactor */
- tcbu->a = BHAMA(fr->nbfp, fr->ntype, ati, atj);
- tcbu->b = BHAMB(fr->nbfp, fr->ntype, ati, atj);
- tcbu->c = BHAMC(fr->nbfp, fr->ntype, ati, atj)/6.0;
- }
-
- for (i = 0; (i < tcr->nQ); i++)
- {
- tcq = &(tcr->tcQ[i]);
- for (j = 0; (j < md->nr); j++)
- {
- if (md->typeA[j] == tcq->at_i)
- {
- tcr->tcQ[i].Q = md->chargeA[j];
- break;
- }
- }
- }
- for (i = 0; (i < tcr->nIP); i++)
- {
- /* Let's just copy the whole struct !*/
- type = tcr->tIP[i].type;
- tcr->tIP[i].iprint = idef->iparams[type];
- }
-}
-
-void write_gct(const char *fn, t_coupl_rec *tcr, t_idef *idef)
-{
- FILE *fp;
- int i, ftype;
-
- fp = gmx_fio_fopen(fn, "w");
- nice_header(fp, fn);
- fprintf(fp, "%-15s = %12g ; Reference pressure for coupling\n",
- eoNames[eoPres], tcr->ref_value[eoPres]);
- fprintf(fp, "%-15s = %12g ; Reference potential energy\n",
- eoNames[eoEpot], tcr->ref_value[eoEpot]);
- fprintf(fp, "%-15s = %12g ; Reference distance\n",
- eoNames[eoDist], tcr->ref_value[eoDist]);
- fprintf(fp, "%-15s = %12g ; Reference dipole\n",
- eoNames[eoMu], tcr->ref_value[eoMu]);
- fprintf(fp, "%-15s = %12g ; Reference force\n",
- eoNames[eoForce], tcr->ref_value[eoForce]);
- fprintf(fp, "%-15s = %12g ; Reference force in X dir\n",
- eoNames[eoFx], tcr->ref_value[eoFx]);
- fprintf(fp, "%-15s = %12g ; Reference force in Y dir\n",
- eoNames[eoFy], tcr->ref_value[eoFy]);
- fprintf(fp, "%-15s = %12g ; Reference force in Z dir\n",
- eoNames[eoFz], tcr->ref_value[eoFz]);
- fprintf(fp, "%-15s = %12g ; Reference pres in X dir\n",
- eoNames[eoPx], tcr->ref_value[eoPx]);
- fprintf(fp, "%-15s = %12g ; Reference pres in Y dir\n",
- eoNames[eoPy], tcr->ref_value[eoPy]);
- fprintf(fp, "%-15s = %12g ; Reference pres in Z dir\n",
- eoNames[eoPz], tcr->ref_value[eoPz]);
- fprintf(fp, "%-15s = %12g ; Polarizability used for the Epot correction\n",
- eoNames[eoPolarizability], tcr->ref_value[eoPolarizability]);
- fprintf(fp, "%-15s = %12g ; Gas phase dipole moment used for Epot correction\n",
- eoNames[eoDipole], tcr->ref_value[eoDipole]);
- fprintf(fp, "%-15s = %12d ; Memory for coupling. Makes it converge faster.\n",
- eoNames[eoMemory], tcr->nmemory);
- fprintf(fp, "%-15s = %12s ; Use intermolecular Epot only (LJ+Coul)\n",
- eoNames[eoInter], yesno_names[tcr->bInter]);
- fprintf(fp, "%-15s = %12s ; Use virial iso pressure\n",
- eoNames[eoUseVirial], yesno_names[tcr->bVirial]);
- fprintf(fp, "%-15s = %12d ; Combination rule, same coding as in grompp.\n",
- eoNames[eoCombRule], tcr->combrule);
-
- fprintf(fp, "\n; Q-Coupling %6s %12s\n", "type", "xi");
- for (i = 0; (i < tcr->nQ); i++)
- {
- fprintf(fp, "%-8s = %8s %6d %12g\n",
- "Q", eoNames[tcr->tcQ[i].eObs], tcr->tcQ[i].at_i, tcr->tcQ[i].xi_Q);
- }
-
- fprintf(fp, "\n; %8s %8s %6s %6s %12s %12s\n", "Couple", "To",
- "i-type", "j-type", "xi-c6", "xi-c12");
- fprintf(fp, "; j-type == -1 means mixing rules will be applied!\n");
- for (i = 0; (i < tcr->nLJ); i++)
- {
- fprintf(fp, "%-8s = %8s %6d %6d %12g %12g\n",
- "LJ", eoNames[tcr->tcLJ[i].eObs],
- tcr->tcLJ[i].at_i, tcr->tcLJ[i].at_j,
- tcr->tcLJ[i].xi_6, tcr->tcLJ[i].xi_12);
- }
-
- fprintf(fp, "\n; %8s %8s %6s %6s %12s %12s %12s\n", "Couple", "To",
- "i-type", "j-type", "xi-A", "xi-B", "xi-C");
- fprintf(fp, "; j-type == -1 means mixing rules will be applied!\n");
- for (i = 0; (i < tcr->nBU); i++)
- {
- fprintf(fp, "%-8s = %8s %6d %6d %12g %12g %12g\n",
- "BU", eoNames[tcr->tcBU[i].eObs],
- tcr->tcBU[i].at_i, tcr->tcBU[i].at_j,
- tcr->tcBU[i].xi_a, tcr->tcBU[i].xi_b, tcr->tcBU[i].xi_c);
- }
-
- fprintf(fp, "\n; More Coupling\n");
- for (i = 0; (i < tcr->nIP); i++)
- {
- ftype = idef->functype[tcr->tIP[i].type];
- switch (ftype)
- {
- case F_BONDS:
- fprintf(fp, "%-15s = %-8s %4d %12g %12g\n",
- "Bonds", eoNames[tcr->tIP[i].eObs], tcr->tIP[i].type,
- tcr->tIP[i].xi.harmonic.krA,
- tcr->tIP[i].xi.harmonic.rA);
- break;
- default:
- fprintf(stderr, "ftype %s not supported (yet)\n",
- interaction_function[ftype].longname);
- }
- }
- gmx_fio_fclose(fp);
-}
-
-static gmx_bool add_lj(int *nLJ, t_coupl_LJ **tcLJ, char *s, gmx_bool bObsUsed[])
-{
- int j, ati, atj, eo;
- char buf[256];
- double xi6, xi12;
-
- if (sscanf(s, "%s%d%d%lf%lf", buf, &ati, &atj, &xi6, &xi12) != 5)
- {
- return TRUE;
- }
- if ((eo = Name2eo(buf)) == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for LJ coupling: %s", buf);
- }
-
- for (j = 0; (j < *nLJ); j++)
- {
- if ((((*tcLJ)[j].at_i == ati) && ((*tcLJ)[j].at_j == atj)) &&
- ((*tcLJ)[j].xi_6 || (*tcLJ)[j].xi_12) &&
- ((*tcLJ)[j].eObs == eo))
- {
- break;
- }
- }
- if (j == *nLJ)
- {
- ++(*nLJ);
- srenew((*tcLJ), *nLJ);
- }
- else
- {
- fprintf(stderr, "\n*** WARNING: overwriting entry for LJ coupling '%s'\n", s);
- }
-
- clear_lj(&((*tcLJ)[j]));
- if (((*tcLJ)[j].eObs = eo) == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for LJ coupling: %s", buf);
- }
- (*tcLJ)[j].at_i = ati;
- (*tcLJ)[j].at_j = atj;
- (*tcLJ)[j].xi_6 = xi6;
- (*tcLJ)[j].xi_12 = xi12;
- bObsUsed[eo] = TRUE;
-
- return FALSE;
-}
-
-static gmx_bool add_bu(int *nBU, t_coupl_BU **tcBU, char *s, gmx_bool bObsUsed[])
-{
- int j, ati, atj, eo;
- char buf[256];
- double xia, xib, xic;
-
- if (sscanf(s, "%s%d%d%lf%lf%lf", buf, &ati, &atj, &xia, &xib, &xic) != 6)
- {
- return TRUE;
- }
- if ((eo = Name2eo(buf)) == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for BU coupling: %s", buf);
- }
-
- for (j = 0; (j < *nBU); j++)
- {
- if ((((*tcBU)[j].at_i == ati) && ((*tcBU)[j].at_j == atj)) &&
- ((*tcBU)[j].xi_a || (*tcBU)[j].xi_b || (*tcBU)[j].xi_c ) &&
- ((*tcBU)[j].eObs == eo))
- {
- break;
- }
- }
- if (j == *nBU)
- {
- ++(*nBU);
- srenew((*tcBU), *nBU);
- }
- else
- {
- fprintf(stderr, "\n*** WARNING: overwriting entry for BU coupling '%s'\n", s);
- }
-
- clear_bu(&((*tcBU)[j]));
- if (((*tcBU)[j].eObs = eo) == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for BU coupling: %s", buf);
- }
- (*tcBU)[j].at_i = ati;
- (*tcBU)[j].at_j = atj;
- (*tcBU)[j].xi_a = xia;
- (*tcBU)[j].xi_b = xib;
- (*tcBU)[j].xi_c = xic;
- bObsUsed[eo] = TRUE;
-
- return FALSE;
-}
-
-static gmx_bool add_ip(int *nIP, t_coupl_iparams **tIP, char *s, int ftype, gmx_bool bObsUsed[])
-{
- int i, eo, type;
- char buf[256];
- double kb, b0;
-
- switch (ftype)
- {
- case F_BONDS:
- /* Pick out the type */
- if (sscanf(s, "%s%d", buf, &type) != 2)
- {
- return TRUE;
- }
- if ((eo = Name2eo(buf)) == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for IP coupling: %s", buf);
- }
-
- /* Check whether this entry is there already */
- for (i = 0; (i < *nIP); i++)
- {
- if ((*tIP)[i].type == type)
- {
- break;
- }
- }
- if (i < *nIP)
- {
- fprintf(stderr, "*** WARNING: overwriting entry for type %d\n", type);
- }
- else
- {
- i = *nIP;
- srenew((*tIP), i+1);
- (*nIP)++;
- }
- if (sscanf(s, "%s%d%lf%lf", buf, &type, &kb, &b0) != 4)
- {
- return TRUE;
- }
- (*tIP)[i].type = type;
- (*tIP)[i].eObs = eo;
- (*tIP)[i].xi.harmonic.krA = kb;
- (*tIP)[i].xi.harmonic.rA = b0;
- bObsUsed[eo] = TRUE;
- break;
- default:
- fprintf(stderr, "ftype %s not supported (yet)\n",
- interaction_function[ftype].longname);
- return TRUE;
- }
- return FALSE;
-}
-
-static gmx_bool add_q(int *nQ, t_coupl_Q **tcQ, char *s, gmx_bool bObsUsed[])
-{
- int j, ati, eo;
- char buf[256];
- double xiQ;
-
- if (sscanf(s, "%s%d%lf", buf, &ati, &xiQ) != 3)
- {
- return TRUE;
- }
-
- for (j = 0; (j < *nQ); j++)
- {
- if ((*tcQ)[j].at_i == ati)
- {
- break;
- }
- }
- if (j == *nQ)
- {
- ++(*nQ);
- srenew((*tcQ), *nQ);
- }
- else
- {
- fprintf(stderr, "\n*** WARNING: overwriting entry for Q coupling '%s'\n", s);
- }
-
- clear_q(&((*tcQ)[j]));
- eo = (*tcQ)[j].eObs = Name2eo(buf);
- if ((*tcQ)[j].eObs == -1)
- {
- gmx_fatal(FARGS, "Invalid observable for Q coupling: %s", buf);
- }
- (*tcQ)[j].at_i = ati;
- (*tcQ)[j].xi_Q = xiQ;
- bObsUsed[eo] = TRUE;
-
- return FALSE;
-}
-
-void read_gct(const char *fn, t_coupl_rec *tcr)
-{
- warninp_t wi;
- t_inpfile *inp;
- int i, j, ninp, nQ, nLJ, nBU, nIP;
- gmx_bool bWrong;
-
- wi = init_warning(FALSE, 0);
-
- inp = read_inpfile(fn, &ninp, wi);
-
- for (i = 0; (i < eoObsNR); i++)
- {
- tcr->bObsUsed[i] = FALSE;
- RTYPE (eoNames[i], tcr->ref_value[i], 0.0);
- }
- ITYPE (eoNames[eoMemory], tcr->nmemory, 1);
- ETYPE (eoNames[eoInter], tcr->bInter, yesno_names);
- ETYPE (eoNames[eoUseVirial], tcr->bVirial, yesno_names);
- ITYPE (eoNames[eoCombRule], tcr->combrule, 1);
- tcr->tcLJ = NULL;
- tcr->tcBU = NULL;
- tcr->tcQ = NULL;
- tcr->tIP = NULL;
- nQ = nLJ = nBU = nIP = 0;
-
- for (i = 0; (i < ninp); i++)
- {
- bWrong = FALSE;
- if (gmx_strcasecmp(inp[i].name, "LJ") == 0)
- {
- bWrong = add_lj(&nLJ, &(tcr->tcLJ), inp[i].value, tcr->bObsUsed);
- }
- else if (gmx_strcasecmp(inp[i].name, "BU") == 0)
- {
- bWrong = add_bu(&nBU, &(tcr->tcBU), inp[i].value, tcr->bObsUsed);
- }
- else if (gmx_strcasecmp(inp[i].name, "Q") == 0)
- {
- bWrong = add_q(&nQ, &(tcr->tcQ), inp[i].value, tcr->bObsUsed);
- }
- else if (gmx_strcasecmp(inp[i].name, "Bonds") == 0)
- {
- bWrong = add_ip(&nIP, &(tcr->tIP), inp[i].value, F_BONDS, tcr->bObsUsed);
- }
-
- if (bWrong)
- {
- fprintf(stderr, "Wrong line in %s: '%s = %s'\n",
- fn, inp[i].name, inp[i].value);
- }
- /*sfree(inp[i].name);
- sfree(inp[i].value);*/
- }
- /* Check which ones have to be printed */
- for (i = 1; (i < nQ); i++)
- {
- for (j = 0; (j < i); j++)
- {
- if (tcr->tcQ[i].at_i == tcr->tcQ[j].at_i)
- {
- tcr->tcQ[j].bPrint = FALSE;
- }
- }
- }
- for (i = 1; (i < nLJ); i++)
- {
- for (j = 0; (j < i); j++)
- {
- if (((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_i) &&
- (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_j)) ||
- ((tcr->tcLJ[i].at_i == tcr->tcLJ[j].at_j) &&
- (tcr->tcLJ[i].at_j == tcr->tcLJ[j].at_i)))
- {
- tcr->tcLJ[j].bPrint = FALSE;
- }
- }
- }
-
- for (i = 1; (i < nBU); i++)
- {
- for (j = 0; (j < i); j++)
- {
- if (((tcr->tcBU[i].at_i == tcr->tcBU[j].at_i) &&
- (tcr->tcBU[i].at_j == tcr->tcBU[j].at_j)) ||
- ((tcr->tcBU[i].at_i == tcr->tcBU[j].at_j) &&
- (tcr->tcBU[i].at_j == tcr->tcBU[j].at_i)))
- {
- tcr->tcBU[j].bPrint = FALSE;
- }
- }
- }
-
- tcr->nQ = nQ;
- tcr->nLJ = nLJ;
- tcr->nBU = nBU;
- tcr->nIP = nIP;
-
- sfree(inp);
-
- done_warning(wi, FARGS);
-}
#include "xvgr.h"
#include "physics.h"
#include "names.h"
-#include "xmdrun.h"
+#include "force.h"
#include "ionize.h"
#include "disre.h"
#include "orires.h"
real timestep = 0;
double tcount = 0;
gmx_bool bIonize = FALSE;
- gmx_bool bTCR = FALSE, bConverged = TRUE, bOK, bSumEkinhOld, bExchanged;
+ gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bExchanged;
gmx_bool bAppend;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
int a0, a1, gnx = 0, ii;
atom_id *grpindex = NULL;
char *grpname;
- t_coupl_rec *tcr = NULL;
rvec *xcopy = NULL, *vcopy = NULL, *cbuf = NULL;
matrix boxcopy = {{0}}, lastbox;
tensor tmpvir;
}
}
- /* Check whether we have to GCT stuff */
- bTCR = ftp2bSet(efGCT, nfile, fnm);
- if (bTCR)
- {
- if (MASTER(cr))
- {
- fprintf(stderr, "Will do General Coupling Theory!\n");
- }
- gnx = top_global->mols.nr;
- snew(grpindex, gnx);
- for (i = 0; (i < gnx); i++)
- {
- grpindex[i] = i;
- }
- }
-
if (repl_ex_nst > 0)
{
/* We need to be sure replica exchange can only occur
(bNS ? GMX_FORCE_NS : 0) | force_flags);
}
- if (bTCR)
- {
- mu_aver = calc_mu_aver(cr, state->x, mdatoms->chargeA,
- mu_tot, &top_global->mols, mdatoms, gnx, grpindex);
- }
-
- if (bTCR && bFirstStep)
- {
- tcr = init_coupling(fplog, nfile, fnm, cr, fr, mdatoms, &(top->idef));
- fprintf(fplog, "Done init_coupling\n");
- fflush(fplog);
- }
-
if (bVV && !bStartingFromCpt && !bRerunMD)
/* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
{
bSumEkinhOld = TRUE;
}
- if (bTCR)
- {
- /* Only do GCT when the relaxation of shells (minimization) has converged,
- * otherwise we might be coupling to bogus energies.
- * In parallel we must always do this, because the other sims might
- * update the FF.
- */
-
- /* Since this is called with the new coordinates state->x, I assume
- * we want the new box state->box too. / EL 20040121
- */
- do_coupling(fplog, oenv, nfile, fnm, tcr, t, step, enerd->term, fr,
- ir, MASTER(cr),
- mdatoms, &(top->idef), mu_aver,
- top_global->mols.nr, cr,
- state->box, total_vir, pres,
- mu_tot, state->x, f, bConverged);
- debug_gmx();
- }
-
/* ######### BEGIN PREPARING EDR OUTPUT ########### */
/* use the directly determined last velocity, not actually the averaged half steps */
{ efXVG, "-tpid", "tpidist", ffOPTWR },
{ efEDI, "-ei", "sam", ffOPTRD },
{ efXVG, "-eo", "edsam", ffOPTWR },
- { efGCT, "-j", "wham", ffOPTRD },
- { efGCT, "-jo", "bam", ffOPTWR },
- { efXVG, "-ffout", "gct", ffOPTWR },
{ efXVG, "-devout", "deviatie", ffOPTWR },
{ efXVG, "-runav", "runaver", ffOPTWR },
{ efXVG, "-px", "pullx", ffOPTWR },
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
- */
-
-#ifndef _xmdrun_h
-#define _xmdrun_h
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include "typedefs.h"
-#include "network.h"
-#include "tgroup.h"
-#include "filenm.h"
-#include "mshift.h"
-#include "force.h"
-#include <time.h>
-#include "edsam.h"
-#include "mdebin.h"
-#include "vcm.h"
-#include "vsite.h"
-#include "gmx_wallcycle.h"
-
-/* This file contains XMDRUN datatypes and function prototypes, grouped
- * neatly according to parts of the functionalisty
- */
-
-/* GENERAL COUPLING THEORY (GCT) STUFF */
-enum {
- eoPres, eoEpot, eoVir, eoDist, eoMu, eoForce, eoFx, eoFy, eoFz,
- eoPx, eoPy, eoPz,
- eoPolarizability, eoDipole, eoObsNR,
- eoMemory = eoObsNR, eoInter, eoUseVirial, eoCombRule, eoNR
-};
-extern const char *eoNames[eoNR];
-
-typedef struct {
- int at_i, at_j; /* Atom type # for i and j */
- int eObs; /* Observable to couple to */
- gmx_bool bPrint; /* Does this struct have to be printed */
- real c6, c12; /* Actual value of params */
- real xi_6, xi_12; /* Constants for coupling C6 and C12 */
-} t_coupl_LJ;
-
-typedef struct {
- int at_i, at_j; /* Atom type # for i and j */
- int eObs; /* Observable to couple to */
- gmx_bool bPrint; /* Does this struct have to be printed */
- real a, b, c; /* Actual value of params */
- real xi_a, xi_b, xi_c; /* Constants for coupling A, B and C */
-} t_coupl_BU;
-
-typedef struct {
- int at_i; /* Atom type */
- int eObs; /* Observable to couple to */
- gmx_bool bPrint; /* Does this struct have to be printed */
- real Q; /* Actual value of charge */
- real xi_Q; /* Constant for coupling Q */
-} t_coupl_Q;
-
-typedef struct {
- int type; /* Type number in the iparams struct */
- int eObs; /* Observable to couple to */
- t_iparams xi; /* Parameters that need to be changed */
- t_iparams iprint;
-} t_coupl_iparams;
-
-typedef struct {
- real act_value[eoObsNR];
- real av_value [eoObsNR];
- real ref_value[eoObsNR];
- gmx_bool bObsUsed[eoObsNR];
- int nLJ, nBU, nQ, nIP;
- t_coupl_LJ *tcLJ;
- t_coupl_BU *tcBU;
- t_coupl_Q *tcQ;
- t_coupl_iparams *tIP;
- int nmemory;
- gmx_bool bInter;
- gmx_bool bVirial;
- int combrule;
-} t_coupl_rec;
-
-extern void write_gct(const char *fn, t_coupl_rec *tcr, t_idef *idef);
-
-extern void read_gct(const char *fn, t_coupl_rec *tcr);
-
-extern void comm_tcr(FILE *log, t_commrec *cr, t_coupl_rec **tcr);
-
-extern void copy_ff(t_coupl_rec *tcr, t_forcerec *fr, t_mdatoms *md,
- t_idef *idef);
-
-extern t_coupl_rec *init_coupling(FILE *log, int nfile, const t_filenm fnm[],
- t_commrec *cr, t_forcerec *fr, t_mdatoms *md,
- t_idef *idef);
-
-extern void calc_force(int natom, rvec f[], rvec fff[]);
-
-extern void calc_f_dev(int natoms, real charge[], rvec x[], rvec f[],
- t_idef *idef, real *xiH, real *xiS);
-
-extern void do_coupling(FILE *log, const output_env_t oenv, int nfile,
- const t_filenm fnm[],
- t_coupl_rec *tcr, real t, int step, real ener[],
- t_forcerec *fr, t_inputrec *ir, gmx_bool bMaster,
- t_mdatoms *md, t_idef *idef, real mu_aver, int nmols,
- t_commrec *cr, matrix box, tensor virial,
- tensor pres, rvec mu_tot,
- rvec x[], rvec f[], gmx_bool bDoIt);
-
-/* CODE TO ADD SPECIAL 2-DIMENSIONAL LENNARD-JONES CORRECTION TO FORCES AND ENERGY */
-extern void do_glas(FILE *log, int start, int homenr, rvec x[], rvec f[],
- t_forcerec *fr, t_mdatoms *md, int atnr, t_inputrec *ir,
- real ener[]);
-
-extern real mol_dipole(int k0, int k1, rvec x[], real q[]);
-/* Calculate total dipole for group of atoms */
-
-extern real calc_mu_aver(t_commrec *cr, rvec x[], real q[], rvec mu,
- t_block *mols, t_mdatoms *md, int gnx, atom_id grpindex[]);
-/* Compute average dipole */
-
-#endif /* _xmdrun_h */
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
- */
-#include "typedefs.h"
-#include "xmdrun.h"
-#include "vec.h"
-
-real mol_dipole(int k0, int k1, rvec x[], real q[])
-{
- int k, m;
- rvec mu;
-
- clear_rvec(mu);
- for (k = k0; (k < k1); k++)
- {
- for (m = 0; (m < DIM); m++)
- {
- mu[m] += q[k]*x[k][m];
- }
- }
- return norm(mu); /* Dipole moment of this molecule in e nm */
-}
-
-real calc_mu_aver(t_commrec *cr, rvec x[], real q[], rvec mu,
- t_block *mols, t_mdatoms *md, int gnx, atom_id grpindex[])
-{
- int i, start, end;
- real mu_ave;
-
- start = md->start;
- end = md->homenr + start;
-
- /*
- clear_rvec(mu);
- for(i=start; (i<end); i++)
- for(m=0; (m<DIM); m++)
- mu[m] += q[i]*x[i][m];
- if (PAR(cr)) {
- gmx_sum(DIM,mu,cr);
- }
- */
- /* I guess we have to parallelise this one! */
-
- if (gnx > 0)
- {
- mu_ave = 0.0;
- for (i = 0; (i < gnx); i++)
- {
- int gi = grpindex[i];
- mu_ave += mol_dipole(mols->index[gi], mols->index[gi+1], x, q);
- }
-
- return(mu_ave/gnx);
- }
- else
- {
- return 0;
- }
-}
-
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff