2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2006, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
55 real ener(matrix P,real e,real e0,int nmol,real kp,real ke,gmx_bool bPScal)
58 return (kp*(sqr(P[XX][XX]+P[YY][YY]+P[ZZ][ZZ]-3))+
61 return (kp*(sqr(P[XX][XX]-1)+sqr(P[YY][YY]-1)+sqr(P[ZZ][ZZ]-1)+
62 sqr(P[XX][YY])+sqr(P[XX][ZZ])+sqr(P[YY][ZZ])) +
66 void do_sim(char *enx,
67 t_topology *top,rvec *x,rvec *v,t_inputrec *ir,matrix box)
69 char *tpx = "optwat.tpr";
72 write_tpx(tpx,0,0.0,0.0,ir,box,top->atoms.nr,x,v,NULL,top);
73 sprintf(buf,"xmdrun -s %s -e %s >& /dev/null",tpx,enx);
77 void get_results(char *enx,real P[],real *epot,int pindex,int eindex)
82 gmx_enxnm_t *nms=NULL;
85 fp_ene = open_enx(enx,"r");
87 do_enxnms(fp_ene,&nre,&nms);
89 /* Read until the last frame */
90 while (do_enx(fp_ene,&fr));
94 *epot = fr.ener[eindex].e;
95 for(i=pindex; (i<pindex+9); i++)
96 P[i-pindex] = fr.ener[i].e;
101 void copy_iparams(int nr,t_iparams dest[],t_iparams src[])
103 memcpy(dest,src,nr*sizeof(dest[0]));
106 void rand_step(FILE *fp,int nr,t_iparams ip[],int *seed,real frac)
112 i = (int) (rando(seed)*nr);
113 } while (ip[i].lj.c12 == 0.0);
116 ff = frac*rando(seed);
119 if (rando(seed) > 0.5) {
120 ip[i].lj.c12 *= 1.0+ff;
121 fprintf(fp,"Increasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
124 ip[i].lj.c12 *= 1.0-ff;
125 fprintf(fp,"Decreasing c12[%d] by %g%% to %g\n",i,100*ff,ip[i].lj.c12);
129 void pr_progress(FILE *fp,int nit,tensor P,real epot,real eFF,
130 double mc_crit,gmx_bool bConv,gmx_bool bAccept)
132 fprintf(fp,"Iter %3d, eFF = %g, Converged = %s, Accepted = %s\n",
133 nit,eFF,yesno_names[bConv],yesno_names[bAccept]);
134 fprintf(fp,"Epot = %g Pscal = %g, mc_crit = %g\n",epot,
136 pr_rvecs(fp,0,"Pres",P,DIM);
137 fprintf(fp,"-----------------------------------------------------\n");
141 int main(int argc,char *argv[])
143 static char *desc[] = {
144 "[TT]optwat[tt] optimizes the force field parameter set of a molecular crystal",
145 "to reproduce the pressure tensor and experimental energy.[PAR]",
146 "Note that for good results the [TT].tpx[tt] file must contain input for a",
147 "simulated annealing run, or a single point energy calculation at 0 K"
150 { efTPX, NULL, NULL, ffREAD },
151 { efEDR, "-e", NULL, ffRW },
152 { efLOG, "-g", NULL, ffWRITE }
154 #define NFILE asize(fnm)
156 static real epot0 = -57, tol = 1, kT = 0.0;
157 static real kp = 1, ke = 100, frac = 0.1;
158 static int maxnit = 100, eindex = 5, pindex = 19;
159 static int seed = 1993;
160 static gmx_bool bPScal = FALSE;
161 static t_pargs pa[] = {
162 { "-epot0", FALSE, etREAL, {&epot0},
163 "Potential energy in kJ/mol" },
164 { "-tol", FALSE, etREAL, {&tol},
165 "Tolerance for converging" },
166 { "-nit", FALSE, etINT, {&maxnit},
167 "Max number of iterations" },
168 { "-seed", FALSE, etINT, {&seed},
169 "Random seed for MC steps" },
170 { "-frac", FALSE, etREAL, {&frac},
171 "Maximum fraction by which to change parameters. Actual fraction is random between 0 and this parameter" },
172 { "-pindex", FALSE, etINT, {&pindex},
173 "Index of P[X][X] in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
174 { "-eindex", FALSE, etINT, {&pindex},
175 "Index of Epot in the energy file (check with [TT]g_energy[tt] and subtract 1)" },
176 { "-kp", FALSE, etREAL, {&kp},
177 "Force constant for pressure components"},
178 { "-ke", FALSE, etREAL, {&ke},
179 "Force constant for energy component"},
180 { "-kT", FALSE, etREAL, {&kT},
181 "Boltzmann Energy for Monte Carlo" },
182 { "-pscal", FALSE, etBOOL, {&bPScal},
183 "Optimize params for scalar pressure, instead of tensor" }
195 int i,step,natoms,nmol,nit,atnr2;
196 real t,lambda,epot,eFF[2];
198 gmx_bool bConverged,bAccept;
201 CopyRight(stdout,argv[0]);
202 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
203 asize(desc),desc,0,NULL);
205 /* Read initial topology and coordaintes etc. */
206 read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&sh,TRUE,NULL,NULL);
209 read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,box,&natoms,
212 /* Open log file and print options */
213 fp = ftp2FILE(efLOG,NFILE,fnm,"w");
214 fprintf(fp,"%s started with the following parameters\n",argv[0]);
215 fprintf(fp,"epot = %8g ke = %8g kp = %8g\n",epot0,ke,kp);
216 fprintf(fp,"maxnit = %8d tol = %8g seed = %8d\n",maxnit,tol,seed);
217 fprintf(fp,"frac = %8g pindex = %8d eindex = %8d\n",frac,pindex,eindex);
218 fprintf(fp,"kT = %8g pscal = %8s\n",kT,gmx_bool_names[bPScal]);
220 /* Unpack some topology numbers */
221 nmol = top.blocks[ebMOLS].nr;
222 atnr2 = top.idef.atnr*top.idef.atnr;
224 /* Copy input params */
226 snew(ip[next],atnr2);
227 copy_iparams(atnr2,ip[cur],top.idef.iparams);
228 copy_iparams(atnr2,ip[next],top.idef.iparams);
230 /* Loop over iterations */
235 rand_step(fp,atnr2,ip[next],&seed,frac);
236 copy_iparams(atnr2,top.idef.iparams,ip[next]);
238 do_sim(ftp2fn(efEDR,NFILE,fnm),&top,xx,vv,&ir,box);
240 get_results(ftp2fn(efEDR,NFILE,fnm),P[0],&epot,pindex,eindex);
242 /* Calculate penalty */
243 eFF[(nit > 0) ? next : cur] = ener(P,epot,epot0,nmol,kp,ke,bPScal);
245 bConverged = (eFF[(nit > 0) ? next : cur] < tol);
248 /* Do Metropolis criterium */
250 mc_crit = exp(-(eFF[next]-eFF[cur])/kT);
251 bAccept = ((eFF[next] < eFF[cur]) ||
252 ((kT > 0) && (mc_crit > rando(&seed))));
253 pr_progress(fp,nit,P,epot/nmol,eFF[next],mc_crit,
260 /* Restore old parameters */
261 copy_iparams(atnr2,ip[next],ip[cur]);
265 pr_progress(fp,nit,P,epot/nmol,eFF[cur],mc_crit,bConverged,FALSE);
268 } while (!bConverged && (nit < maxnit));
270 for(i=0; (i<atnr2); i++)
271 pr_iparams(fp,F_LJ,&ip[cur][i]);