#define MD_IONIZE (1<<3)
#define MD_RERUN (1<<4)
#define MD_RERUN_VSITE (1<<5)
-#define MD_FFSCAN (1<<6)
#define MD_SEPPOT (1<<7)
#define MD_PARTDEC (1<<9)
#define MD_DDBONDCHECK (1<<10)
+++ /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
- */
-/**H*O*C**************************************************************
-** **
-** No.!Version! Date ! Request ! Modification ! Author **
-** ---+-------+------+---------+---------------------------+------- **
-** 1 + 3.1 +5/18/95+ - + strategy DE/rand-to-best/1+ Storn **
-** + + + + included + **
-** 1 + 3.2 +6/06/95+C.Fleiner+ change loops into memcpy + Storn **
-** 2 + 3.2 +6/06/95+ - + update comments + Storn **
-** 1 + 3.3 +6/15/95+ K.Price + strategy DE/best/2 incl. + Storn **
-** 2 + 3.3 +6/16/95+ - + comments and beautifying + Storn **
-** 3 + 3.3 +7/13/95+ - + upper and lower bound for + Storn **
-** + + + + initialization + **
-** 1 + 3.4 +2/12/96+ - + increased printout prec. + Storn **
-** 1 + 3.5 +5/28/96+ - + strategies revisited + Storn **
-** 2 + 3.5 +5/28/96+ - + strategy DE/rand/2 incl. + Storn **
-** 1 + 3.6 +8/06/96+ K.Price + Binomial Crossover added + Storn **
-** 2 + 3.6 +9/30/96+ K.Price + cost variance output + Storn **
-** 3 + 3.6 +9/30/96+ - + alternative to ASSIGND + Storn **
-** 4 + 3.6 +10/1/96+ - + variable checking inserted + Storn **
-** 5 + 3.6 +10/1/96+ - + strategy indic. improved + Storn **
-** **
-***H*O*C*E***********************************************************/
-
-/* Adopted for use in GROMACS by David van der Spoel, Oct 2001 */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <memory.h>
-#include "typedefs.h"
-#include "smalloc.h"
-#include "futil.h"
-#include "genalg.h"
-#include "gmx_fatal.h"
-#include "random.h"
-#include "txtdump.h"
-#include "vec.h"
-#include "main.h"
-#include "gmxfio.h"
-
-static const char *strat[] = {
- "DE/best/1/exp", "DE/rand/1/exp",
- "DE/rand-to-best/1/exp", "DE/best/2/exp",
- "DE/rand/2/exp", "DE/best/1/bin",
- "DE/rand/1/bin", "DE/rand-to-best/1/bin",
- "DE/best/2/bin", "DE/rand/2/bin"
-};
-
-/*------------------------Macros----------------------------------------*/
-static void assignd(int D, real a[], real b[])
-{
- int i;
-
- for (i = 0; (i < D); i++)
- {
- a[i] = b[i];
- }
-}
-
-static real **make2d(int n, int m)
-{
- int i;
- real **r;
-
- snew(r, n);
- for (i = 0; (i < n); i++)
- {
- snew(r[i], m);
- }
- return r;
-}
-
-static void copy2range(int D, real c[], t_range r[])
-{
- int i;
-
- for (i = 0; (i < D); i++)
- {
- /* Range check */
- while (c[i] < r[i].rmin)
- {
- c[i] += r[i].dr;
- }
- while (c[i] > r[i].rmax)
- {
- c[i] -= r[i].dr;
- }
- /* if (c[i] < r[i].rmin)
- c[i] = r[i].rmin;
- if (c[i] > r[i].rmax)
- c[i] = r[i].rmax;
- */
- r[i].rval = c[i];
- }
-}
-
-t_genalg *init_ga(FILE *fplog, const char *infile, int D, t_range range[])
-{
- FILE *fpin_ptr;
- t_genalg *ga;
- double ff = 0.0, cr = 0.0;
- int i, j;
-
- /*------Initializations----------------------------*/
- snew(ga, 1);
- /*-----Read input data------------------------------------------------*/
- fpin_ptr = gmx_fio_fopen(infile, "r");
- if (fscanf(fpin_ptr, "%d", &ga->NP) != 1 || /*---choice of strategy---*/
- fscanf(fpin_ptr, "%d", &ga->strategy) != 1 || /*---choice of strategy---*/
- fscanf(fpin_ptr, "%lf", &ff) != 1 || /*---weight factor------------*/
- fscanf(fpin_ptr, "%lf", &cr) != 1 || /*---crossing over factor-----*/
- fscanf(fpin_ptr, "%d", &ga->seed) != 1 || /*---random seed----------*/
- gmx_fio_fclose(fpin_ptr) != 1)
- {
- gmx_fatal(FARGS, "Error reading from file %s", infile);
- }
-
- ga->FF = ff;
- ga->CR = cr;
- ga->D = D;
- ga->ipop = 0;
- ga->gen = 0;
-
- /* Allocate memory */
- ga->pold = make2d(ga->NP, ga->D);
- ga->pnew = make2d(ga->NP, ga->D);
- snew(ga->tmp, ga->D);
- snew(ga->best, ga->D);
- snew(ga->bestit, ga->D);
- snew(ga->cost, ga->NP);
- snew(ga->msf, ga->NP);
- snew(ga->pres, ga->NP);
- snew(ga->scale, ga->NP);
- snew(ga->energy, ga->NP);
-
- /*-----Checking input variables for proper range--------------*/
- if ((ga->CR < 0) || (ga->CR > 1.0))
- {
- gmx_fatal(FARGS, "CR=%f, should be ex [0,1]", ga->CR);
- }
- if (ga->seed <= 0)
- {
- gmx_fatal(FARGS, "seed=%d, should be > 0", ga->seed);
- }
- if ((ga->strategy < 0) || (ga->strategy > 10))
- {
- gmx_fatal(FARGS, "strategy=%d, should be ex {1-10}", ga->strategy);
- }
-
- /* spread initial population members */
- for (i = 0; (i < ga->NP); i++)
- {
- for (j = 0; (j < ga->D); j++)
- {
- ga->pold[i][j] = value_rand(&(range[j]), &ga->seed);
- }
- }
-
- fprintf(fplog, "-----------------------------------------------\n");
- fprintf(fplog, "Genetic algorithm parameters\n");
- fprintf(fplog, "-----------------------------------------------\n");
- fprintf(fplog, "Number of variables: %d\n", ga->D);
- fprintf(fplog, "Population size: %d\n", ga->NP);
- fprintf(fplog, "Strategy: %s\n", strat[ga->strategy]);
- fprintf(fplog, "Weight factor: %g\n", ga->FF);
- fprintf(fplog, "Crossing over factor: %g\n", ga->CR);
- fprintf(fplog, "Random seed: %d\n", ga->seed);
- fprintf(fplog, "-----------------------------------------------\n");
-
- return ga;
-}
-
-void update_ga(FILE *fpout_ptr, t_range range[], t_genalg *ga)
-{
- static int i_init = 0; /* Initialisation related stuff */
- int i, j, L, n; /* counting variables */
- int r1, r2, r3, r4, r5; /* placeholders for random indexes */
-
- if (i_init < ga->NP)
- {
- /* Copy data for first force evaluation to range array */
- copy2range(ga->D, ga->pold[i_init], range);
-
- i_init++;
- return;
- }
- else
- {
- /* Now starts real genetic stuff, a new trial set is made */
- if (ga->ipop == ga->NP)
- {
- ga->gen++;
- i = ga->ipop = 0;
- }
- else
- {
- i = ga->ipop;
- }
-
- do /* Pick a random population member */
- { /* Endless loop for ga->NP < 2 !!! */
- r1 = (int)(rando(&ga->seed)*ga->NP);
- }
- while (r1 == i);
-
- do /* Pick a random population member */
- { /* Endless loop for ga->NP < 3 !!! */
- r2 = (int)(rando(&ga->seed)*ga->NP);
- }
- while ((r2 == i) || (r2 == r1));
-
- do
- {
- /* Pick a random population member */
- /* Endless loop for ga->NP < 4 !!! */
- r3 = (int)(rando(&ga->seed)*ga->NP);
- }
- while ((r3 == i) || (r3 == r1) || (r3 == r2));
-
- do
- {
- /* Pick a random population member */
- /* Endless loop for ga->NP < 5 !!! */
- r4 = (int)(rando(&ga->seed)*ga->NP);
- }
- while ((r4 == i) || (r4 == r1) || (r4 == r2) || (r4 == r3));
-
- do
- {
- /* Pick a random population member */
- /* Endless loop for ga->NP < 6 !!! */
- r5 = (int)(rando(&ga->seed)*ga->NP);
- }
- while ((r5 == i) || (r5 == r1) || (r5 == r2) || (r5 == r3) || (r5 == r4));
-
-
- /* Choice of strategy
- * We have tried to come up with a sensible naming-convention: DE/x/y/z
- * DE : stands for Differential Evolution
- * x : a string which denotes the vector to be perturbed
- * y : number of difference vectors taken for perturbation of x
- * z : crossover method (exp = exponential, bin = binomial)
- *
- * There are some simple rules which are worth following:
- * 1) ga->FF is usually between 0.5 and 1 (in rare cases > 1)
- * 2) ga->CR is between 0 and 1 with 0., 0.3, 0.7 and 1. being worth to
- * be tried first
- * 3) To start off ga->NP = 10*ga->D is a reasonable choice. Increase ga->NP if
- * misconvergence happens.
- * 4) If you increase ga->NP, ga->FF usually has to be decreased
- * 5) When the DE/ga->best... schemes fail DE/rand... usually works and
- * vice versa
- * EXPONENTIAL ga->CROSSOVER
- *-------DE/ga->best/1/exp-------
- *-------Our oldest strategy but still not bad. However, we have found several
- *-------optimization problems where misconvergence occurs.
- */
- assignd(ga->D, ga->tmp, ga->pold[i]);
- n = (int)(rando(&ga->seed)*ga->D);
- L = 0;
-
- switch (ga->strategy)
- {
- case 1:
- /* strategy DE0 (not in our paper) */
- do
- {
- ga->tmp[n] = ga->bestit[n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
- n = (n+1)%ga->D;
- L++;
- }
- while ((rando(&ga->seed) < ga->CR) && (L < ga->D));
- break;
-
- /* DE/rand/1/exp
- * This is one of my favourite strategies. It works especially
- * well when the ga->bestit[]"-schemes experience misconvergence.
- * Try e.g. ga->FF=0.7 and ga->CR=0.5 * as a first guess.
- */
- case 2:
- /* strategy DE1 in the techreport */
- do
- {
- ga->tmp[n] = ga->pold[r1][n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
- n = (n+1)%ga->D;
- L++;
- }
- while ((rando(&ga->seed) < ga->CR) && (L < ga->D));
- break;
-
- /* DE/rand-to-ga->best/1/exp
- * This strategy seems to be one of the ga->best strategies.
- * Try ga->FF=0.85 and ga->CR=1.
- * If you get misconvergence try to increase ga->NP.
- * If this doesn't help you should play around with all three
- * control variables.
- */
- case 3:
- /* similar to DE2 but generally better */
- do
- {
- ga->tmp[n] = ga->tmp[n] + ga->FF*(ga->bestit[n] - ga->tmp[n]) +
- ga->FF*(ga->pold[r1][n]-ga->pold[r2][n]);
- n = (n+1)%ga->D;
- L++;
- }
- while ((rando(&ga->seed) < ga->CR) && (L < ga->D));
- break;
-
- /* DE/ga->best/2/exp is another powerful strategy worth trying */
- case 4:
- do
- {
- ga->tmp[n] = ga->bestit[n] +
- (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
- n = (n+1)%ga->D;
- L++;
- }
- while ((rando(&ga->seed) < ga->CR) && (L < ga->D));
- break;
-
- /*----DE/rand/2/exp seems to be a robust optimizer for many functions-----*/
- case 5:
- do
- {
- ga->tmp[n] = ga->pold[r5][n] +
- (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
- n = (n+1)%ga->D;
- L++;
- }
- while ((rando(&ga->seed) < ga->CR) && (L < ga->D));
- break;
-
- /*===Essentially same strategies but BINOMIAL ga->CROSSOVER===*/
-
- /*-------DE/ga->best/1/bin------*/
- case 6:
- for (L = 0; L < ga->D; L++)
- {
- /* perform D binomial trials */
- if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1)))
- {
- /* change at least one parameter */
- ga->tmp[n] = ga->bestit[n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
- }
- n = (n+1)%ga->D;
- }
- break;
-
- /*-------DE/rand/1/bin------*/
- case 7:
- for (L = 0; L < ga->D; L++)
- {
- /* perform D binomial trials */
- if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1)))
- {
- /* change at least one parameter */
- ga->tmp[n] = ga->pold[r1][n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
- }
- n = (n+1)%ga->D;
- }
- break;
-
- /*-------DE/rand-to-ga->best/1/bin------*/
- case 8:
- for (L = 0; L < ga->D; L++)
- {
- /* perform ga->D binomial trials */
- if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1)))
- {
- /* change at least one parameter */
- ga->tmp[n] = ga->tmp[n] + ga->FF*(ga->bestit[n] - ga->tmp[n]) +
- ga->FF*(ga->pold[r1][n]-ga->pold[r2][n]);
- }
- n = (n+1)%ga->D;
- }
- break;
-
- /*-------DE/ga->best/2/bin------*/
- case 9:
- for (L = 0; L < ga->D; L++)
- {
- /* perform ga->D binomial trials */
- if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1)))
- {
- /* change at least one parameter */
- ga->tmp[n] = ga->bestit[n] +
- (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
- }
- n = (n+1)%ga->D;
- }
- break;
-
- /*-------DE/rand/2/bin-------*/
- default:
- for (L = 0; L < ga->D; L++)
- {
- /* perform ga->D binomial trials */
- if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1)))
- {
- /* change at least one parameter */
- ga->tmp[n] = ga->pold[r5][n] +
- (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
- }
- n = (n+1)%ga->D;
- }
- break;
- }
-
- /*===Trial mutation now in ga->tmp[]. Test how good this choice really was.==*/
- copy2range(ga->D, ga->tmp, range);
- }
-}
-
-gmx_bool print_ga(FILE *fp, t_genalg *ga, real msf, tensor pres, rvec scale,
- real energy, t_range range[], real tol)
-{
- static int nfeval = 0; /* number of function evaluations */
- static gmx_bool bImproved;
- real trial_cost;
- real cvar; /* computes the cost variance */
- real cmean; /* mean cost */
- int i, j;
- real **pswap;
-
- trial_cost = cost(pres, msf, energy);
- if (nfeval < ga->NP)
- {
- ga->cost[nfeval] = trial_cost;
- ga->msf[nfeval] = msf;
- ga->energy[nfeval] = energy;
- copy_mat(pres, ga->pres[nfeval]);
- copy_rvec(scale, ga->scale[nfeval]);
- if (debug)
- {
- pr_rvec(debug, 0, "scale", scale, DIM, TRUE);
- pr_rvec(debug, 0, "pold ", ga->pold[nfeval]+4, DIM, TRUE);
- }
- nfeval++;
- return FALSE;
- }
- /* When we get here we have done an initial evaluation for all
- * animals in the population
- */
- if (ga->ipop == 0)
- {
- bImproved = FALSE;
- }
-
- /* First iteration after first round of trials */
- if (nfeval == ga->NP)
- {
- /* Evaluate who is ga->best */
- ga->imin = 0;
- for (j = 1; (j < ga->NP); j++)
- {
- if (ga->cost[j] < ga->cost[ga->imin])
- {
- ga->imin = j;
- }
- }
- assignd(ga->D, ga->best, ga->pold[ga->imin]);
- /* save best member ever */
- assignd(ga->D, ga->bestit, ga->pold[ga->imin]);
- /* save best member of generation */
- }
-
- if (trial_cost < ga->cost[ga->ipop])
- {
- if (trial_cost < ga->cost[ga->imin])
- {
- /* Was this a new minimum? */
- /* if so, reset cmin to new low...*/
- ga->imin = ga->ipop;
- assignd(ga->D, ga->best, ga->tmp);
- bImproved = TRUE;
- }
- /* improved objective function value ? */
- ga->cost[ga->ipop] = trial_cost;
-
- ga->msf[ga->ipop] = msf;
- ga->energy[ga->ipop] = energy;
- copy_mat(pres, ga->pres[ga->ipop]);
- copy_rvec(scale, ga->scale[ga->ipop]);
-
- assignd(ga->D, ga->pnew[ga->ipop], ga->tmp);
-
- }
- else
- {
- /* replace target with old value */
- assignd(ga->D, ga->pnew[ga->ipop], ga->pold[ga->ipop]);
- }
- /* #define SCALE_DEBUG*/
-#ifdef SCALE_DEBUG
- if (ga->D > 5)
- {
- rvec dscale;
- rvec_sub(ga->scale[ga->imin], &(ga->best[ga->D-3]), dscale);
- if (norm(dscale) > 0)
- {
- pr_rvec(fp, 0, "scale", scale, DIM, TRUE);
- pr_rvec(fp, 0, "best ", &(ga->best[ga->D-3]), DIM, TRUE);
- fprintf(fp, "imin = %d, ipop = %d, nfeval = %d\n", ga->imin,
- ga->ipop, nfeval);
- gmx_fatal(FARGS, "Scale inconsistency");
- }
- }
-#endif
-
- /* Increase population member count */
- ga->ipop++;
-
- /* End mutation loop through population? */
- if (ga->ipop == ga->NP)
- {
- /* Save ga->best population member of current iteration */
- assignd(ga->D, ga->bestit, ga->best);
-
- /* swap population arrays. New generation becomes old one */
- pswap = ga->pold;
- ga->pold = ga->pnew;
- ga->pnew = pswap;
-
- /*----Compute the energy variance (just for monitoring purposes)-----------*/
- /* compute the mean value first */
- cmean = 0.0;
- for (j = 0; (j < ga->NP); j++)
- {
- cmean += ga->cost[j];
- }
- cmean = cmean/ga->NP;
-
- /* now the variance */
- cvar = 0.0;
- for (j = 0; (j < ga->NP); j++)
- {
- cvar += sqr(ga->cost[j] - cmean);
- }
- cvar = cvar/(ga->NP-1);
-
- /*----Output part----------------------------------------------------------*/
- if (1 || bImproved || (nfeval == ga->NP))
- {
- fprintf(fp, "\nGen: %5d\n Cost: %12.3e <Cost>: %12.3e\n"
- " Ener: %8.4f RMSF: %8.3f Pres: %8.1f %8.1f %8.1f (%8.1f)\n"
- " Box-Scale: %15.10f %15.10f %15.10f\n",
- ga->gen, ga->cost[ga->imin], cmean, ga->energy[ga->imin],
- sqrt(ga->msf[ga->imin]), ga->pres[ga->imin][XX][XX],
- ga->pres[ga->imin][YY][YY], ga->pres[ga->imin][ZZ][ZZ],
- trace(ga->pres[ga->imin])/3.0,
- ga->scale[ga->imin][XX], ga->scale[ga->imin][YY],
- ga->scale[ga->imin][ZZ]);
-
- for (j = 0; (j < ga->D); j++)
- {
- fprintf(fp, "\tbest[%d]=%-15.10f\n", j, ga->best[j]);
- }
- if (ga->cost[ga->imin] < tol)
- {
- for (i = 0; (i < ga->NP); i++)
- {
- fprintf(fp, "Animal: %3d Cost:%8.3f Ener: %8.3f RMSF: %8.3f%s\n",
- i, ga->cost[i], ga->energy[i], sqrt(ga->msf[i]),
- (i == ga->imin) ? " ***" : "");
- for (j = 0; (j < ga->D); j++)
- {
- fprintf(fp, "\tParam[%d][%d]=%15.10g\n", i, j, ga->pold[i][j]);
- }
- }
- return TRUE;
- }
- fflush(fp);
- }
- }
- nfeval++;
-
- return FALSE;
-}
+++ /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 _genalg_h
-#define _genalg_h
-
-typedef struct {
- int np; /* Number of points */
- int atype; /* Atom type */
- int ptype; /* Parameter type */
- real rmin, rmax; /* Minimum and maximum value */
- real dr;
- real rval; /* Current value */
-} t_range;
-
-typedef struct {
- int NP, D;
- int strategy;
- int seed;
- int ipop, gen; /* Current population member and generation */
- int imin; /* Member with lowest energy */
- real CR, FF;
- real **pold, **pnew; /* Old and new populations */
- real *best, *bestit, *cost, *tmp, *msf, *energy;
- tensor *pres;
- rvec *scale;
-} t_genalg;
-
-enum {
- eseSIGMA, eseEPSILON, eseBHAMA, eseBHAMB, eseBHAMC,
- eseCELLX, eseCELLY, eseCELLZ, eseNR
-};
-
-extern real value_rand(t_range *r, int *seed);
-
-extern t_genalg *init_ga(FILE *fplog, const char *infile, int D, t_range range[]);
-
-extern void update_ga(FILE *fpout_ptr, t_range range[], t_genalg *ga);
-
-extern gmx_bool print_ga(FILE *fp, t_genalg *ga, real msf, tensor pres, rvec scale,
- real energy, t_range range[], real tol);
-
-extern real cost(tensor P, real MSF, real energy);
-
-#endif
t_graph *graph = NULL;
globsig_t gs;
gmx_rng_t mcrng = NULL;
- gmx_bool bFFscan;
gmx_groups_t *groups;
gmx_ekindata_t *ekind, *ekind_save;
gmx_shellfc_t shellfc;
/* Check for special mdrun options */
bRerunMD = (Flags & MD_RERUN);
bIonize = (Flags & MD_IONIZE);
- bFFscan = (Flags & MD_FFSCAN);
bAppend = (Flags & MD_APPENDFILES);
if (Flags & MD_RESETCOUNTERSHALFWAY)
{
"If you want less energy communication, set nstlist > 3.\n\n");
}
- if (bRerunMD || bFFscan)
+ if (bRerunMD)
{
ir->nstxtcout = 0;
}
{
bufstate = init_bufstate(state);
}
- if (bFFscan)
- {
- snew(xcopy, state->natoms);
- snew(vcopy, state->natoms);
- copy_rvecn(state->x, xcopy, 0, state->natoms);
- copy_rvecn(state->v, vcopy, 0, state->natoms);
- copy_mat(state->box, boxcopy);
- }
/* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
temperature control */
/* Stop Center of Mass motion */
bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
- /* Copy back starting coordinates in case we're doing a forcefield scan */
- if (bFFscan)
- {
- for (ii = 0; (ii < state->natoms); ii++)
- {
- copy_rvec(xcopy[ii], state->x[ii]);
- copy_rvec(vcopy[ii], state->v[ii]);
- }
- copy_mat(boxcopy, state->box);
- }
-
if (bRerunMD)
{
/* for rerun MD always do Neighbour Searching */
}
}
- if (MASTER(cr) && do_log && !bFFscan)
+ if (MASTER(cr) && do_log)
{
print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */
}
mdatoms->start, mdatoms->start+mdatoms->homenr, state->box, cr);
}
- /* Update force field in ffscan program */
- if (bFFscan)
- {
- if (update_forcefield(fplog,
- nfile, fnm, fr,
- mdatoms->nr, state->x, state->box))
- {
- gmx_finalize_par();
-
- exit(0);
- }
- }
-
/* We write a checkpoint at this MD step when:
* either at an NS step when we signalled through gs,
* or at the last step (but not when we do not want confout),
if (shellfc)
{
/* Now is the time to relax the shells */
- count = relax_shell_flexcon(fplog, cr, bVerbose, bFFscan ? step+1 : step,
+ count = relax_shell_flexcon(fplog, cr, bVerbose, step,
ir, bNS, force_flags,
top,
constr, enerd, fcd,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir, vetanew);
- if (!bOK && !bFFscan)
+ if (!bOK)
{
gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
}
debug_gmx();
if (bLastStep && step_rel == ir->nsteps &&
(Flags & MD_CONFOUT) && MASTER(cr) &&
- !bRerunMD && !bFFscan)
+ !bRerunMD)
{
/* x and v have been collected in write_traj,
* because a checkpoint file will always be written
FALSE, bCalcVir,
state->veta);
}
- if (!bOK && !bFFscan)
+ if (!bOK)
{
gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
}
/* #### We now have r(t+dt) and v(t+dt/2) ############# */
/* The coordinates (x) were unshifted in update */
- if (bFFscan && (shellfc == NULL || bConverged))
- {
- if (print_forcefield(fplog, enerd->term, mdatoms->homenr,
- f, NULL, xcopy,
- &(top_global->mols), mdatoms->massT, pres))
- {
- gmx_finalize_par();
-
- fprintf(stderr, "\n");
- exit(0);
- }
- }
if (!bGStat)
{
/* We will not sum ekinh_old,
t_block *mols, t_mdatoms *md, int gnx, atom_id grpindex[]);
/* Compute average dipole */
-/********************************************************************/
-/* Force field scanning stuff */
-typedef struct {
- real tol, f_max, npow, epot, fac_epot, fac_pres, fac_msf, pres;
- int molsize, nmol;
- gmx_bool bComb, bVerbose, bLogEps;
-} t_ffscan;
-
-
-extern gmx_bool update_forcefield(FILE *fplog,
- int nfile, const t_filenm fnm[], t_forcerec *fr,
- int natoms, rvec x[], matrix box);
-/* Modify the parameters. Return TRUE when the scan is finished. */
-
-extern gmx_bool print_forcefield(FILE *fp, real ener[], int natoms, rvec f[],
- rvec fshake[], rvec x[], t_block *mols, real mass[],
- tensor pres);
-/* Print results. Return TRUE when the scan is finished. */
-
-extern void set_ffvars(t_ffscan *ff);
-/* Set variables for force scanning */
-
#endif /* _xmdrun_h */
* And Hey:
* Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <assert.h>
-
#include "typedefs.h"
-#include "smalloc.h"
-#include "strdb.h"
-#include "string2.h"
#include "xmdrun.h"
#include "vec.h"
-#include "genalg.h"
-#include "random.h"
-#include "macros.h"
real mol_dipole(int k0, int k1, rvec x[], real q[])
{
}
}
-/* Lots of global variables! Yummy... */
-static t_ffscan ff;
-
-void set_ffvars(t_ffscan *fff)
-{
- ff = *fff;
-}
-
-real cost(tensor P, real MSF, real E)
-{
- return (ff.fac_msf*MSF+ff.fac_epot*sqr(E-ff.epot)+ff.fac_pres*
- (sqr(P[XX][XX]-ff.pres)+sqr(P[YY][YY]-ff.pres)+sqr(P[ZZ][ZZ]-ff.pres)));
-
-}
-
-static const char *esenm[eseNR] = { "SIG", "EPS", "BHAMA", "BHAMB", "BHAMC", "CELlX", "CELLY", "CELLZ" };
-static int nparm = 0, *param_val = NULL;
-static t_range *range = NULL;
-static t_genalg *ga = NULL;
-static rvec scale = { 1, 1, 1 };
-
-static void init_range(t_range *r, int np, int atype, int ptype,
- real rmin, real rmax)
-{
- if (rmin > rmax)
- {
- gmx_fatal(FARGS, "rmin (%f) > rmax (%f)", rmin, rmax);
- }
- if (np <= 0)
- {
- gmx_fatal(FARGS, "np (%d) should be > 0", np);
- }
- if ((rmax > rmin) && (np <= 1))
- {
- gmx_fatal(FARGS, "If rmax > rmin, np should be > 1");
- }
- if ((ptype < 0) || (ptype >= eseNR))
- {
- gmx_fatal(FARGS, "ptype (%d) should be < %d", ptype, eseNR);
- }
- r->np = np;
- r->atype = atype;
- r->ptype = ptype;
- r->rmin = rmin;
- r->rmax = rmax;
- r->rval = rmin;
- r->dr = r->rmax - r->rmin;
-}
-
-static t_range *read_range(const char *db, int *nrange)
-{
- int nlines, nr, np, i;
- char **lines;
- t_range *range;
- int atype, ptype;
- double rmin, rmax;
-
- nlines = get_file(db, &lines);
- snew(range, nlines);
-
- nr = 0;
- for (i = 0; (i < nlines); i++)
- {
- strip_comment(lines[i]);
- if (sscanf(lines[i], "%d%d%d%lf%lf", &np, &atype, &ptype, &rmin, &rmax) == 5)
- {
- if (ff.bLogEps && (ptype == eseEPSILON) && (rmin <= 0))
- {
- gmx_fatal(FARGS, "When using logarithmic epsilon increments the minimum"
- "value must be > 0");
- }
- init_range(&range[nr], np, atype, ptype, rmin, rmax);
- nr++;
- }
- }
- fprintf(stderr, "found %d variables to iterate over\n", nr);
-
- *nrange = nr;
-
- for (nr = 0; (nr < nlines); nr++)
- {
- sfree(lines[nr]);
- }
- sfree(lines);
-
- return range;
-}
-
-static real value_range(t_range *r, int n)
-{
- real logrmin, logrmax;
-
- if ((n < 0) || (n > r->np))
- {
- gmx_fatal(FARGS, "Value (%d) out of range for value_range (max %d)", n, r->np);
- }
-
- if (r->np == 1)
- {
- r->rval = r->rmin;
- }
- else
- {
- if ((r->ptype == eseEPSILON) && ff.bLogEps)
- {
- logrmin = log(r->rmin);
- logrmax = log(r->rmax);
- r->rval = exp(logrmin + (n*(logrmax-logrmin))/(r->np-1));
- }
- else
- {
- r->rval = r->rmin+(n*(r->dr))/(r->np-1);
- }
- }
- return r->rval;
-}
-
-real value_rand(t_range *r, int *seed)
-{
- real logrmin, logrmax;
- real mr;
-
- if (r->np == 1)
- {
- r->rval = r->rmin;
- }
- else
- {
- mr = rando(seed);
- if ((r->ptype == eseEPSILON) && ff.bLogEps)
- {
- logrmin = log(r->rmin);
- logrmax = log(r->rmax);
- r->rval = exp(logrmin + mr*(logrmax-logrmin));
- }
- else
- {
- r->rval = r->rmin + mr*(r->rmax-r->rmin);
- }
- }
- if (debug)
- {
- fprintf(debug, "type: %s, value: %g\n", esenm[r->ptype], r->rval);
- }
- return r->rval;
-}
-
-static void update_ff(t_forcerec *fr, int nparm, t_range range[], int param_val[])
-{
- static double *sigma = NULL, *eps = NULL, *c6 = NULL, *cn = NULL, *bhama = NULL, *bhamb = NULL, *bhamc = NULL;
- real val, *nbfp;
- int i, j, atnr;
-
- atnr = fr->ntype;
- nbfp = fr->nbfp;
-
- if (fr->bBHAM)
- {
- if (bhama == NULL)
- {
- assert(bhamb == NULL && bhamc == NULL);
- snew(bhama, atnr);
- snew(bhamb, atnr);
- snew(bhamc, atnr);
- }
- }
- else
- {
- if (sigma == NULL)
- {
- assert(eps == NULL && c6 == NULL && cn == NULL);
- snew(sigma, atnr);
- snew(eps, atnr);
- snew(c6, atnr);
- snew(cn, atnr);
- }
- }
- /* Get current values for everything */
- for (i = 0; (i < nparm); i++)
- {
- if (ga)
- {
- val = range[i].rval;
- }
- else
- {
- val = value_range(&range[i], param_val[i]);
- }
- if (debug)
- {
- fprintf(debug, "val = %g\n", val);
- }
- switch (range[i].ptype)
- {
- case eseSIGMA:
- assert(!fr->bBHAM);
- sigma[range[i].atype] = val;
- break;
- case eseEPSILON:
- assert(!fr->bBHAM);
- eps[range[i].atype] = val;
- break;
- case eseBHAMA:
- assert(fr->bBHAM);
- bhama[range[i].atype] = val;
- break;
- case eseBHAMB:
- assert(fr->bBHAM);
- bhamb[range[i].atype] = val;
- break;
- case eseBHAMC:
- assert(fr->bBHAM);
- bhamc[range[i].atype] = val;
- break;
- case eseCELLX:
- scale[XX] = val;
- break;
- case eseCELLY:
- scale[YY] = val;
- break;
- case eseCELLZ:
- scale[ZZ] = val;
- break;
- default:
- gmx_fatal(FARGS, "Unknown ptype");
- }
- }
- if (fr->bBHAM)
- {
- for (i = 0; (i < atnr); i++)
- {
- for (j = 0; (j <= i); j++)
- {
- BHAMA(nbfp, atnr, i, j) = BHAMA(nbfp, atnr, j, i) = sqrt(bhama[i]*bhama[j]);
- BHAMB(nbfp, atnr, i, j) = BHAMB(nbfp, atnr, j, i) = sqrt(bhamb[i]*bhamb[j]);
- BHAMC(nbfp, atnr, i, j) = BHAMC(nbfp, atnr, j, i) = sqrt(bhamc[i]*bhamc[j]);
- }
- }
- }
- else
- {
- /* Now build a new matrix */
- for (i = 0; (i < atnr); i++)
- {
- c6[i] = 4*eps[i]*pow(sigma[i], 6.0);
- cn[i] = 4*eps[i]*pow(sigma[i], ff.npow);
- }
- for (i = 0; (i < atnr); i++)
- {
- for (j = 0; (j <= i); j++)
- {
- C6(nbfp, atnr, i, j) = C6(nbfp, atnr, j, i) = sqrt(c6[i]*c6[j]);
- C12(nbfp, atnr, i, j) = C12(nbfp, atnr, j, i) = sqrt(cn[i]*cn[j]);
- }
- }
- }
-
- if (debug)
- {
- if (!fr->bBHAM)
- {
- for (i = 0; (i < atnr); i++)
- {
- fprintf(debug, "atnr = %2d sigma = %8.4f eps = %8.4f\n", i, sigma[i], eps[i]);
- }
- }
- for (i = 0; (i < atnr); i++)
- {
- for (j = 0; (j < atnr); j++)
- {
- if (fr->bBHAM)
- {
- fprintf(debug, "i: %2d j: %2d A: %10.5e B: %10.5e C: %10.5e\n", i, j,
- BHAMA(nbfp, atnr, i, j), BHAMB(nbfp, atnr, i, j), BHAMC(nbfp, atnr, i, j));
- }
- else
- {
- fprintf(debug, "i: %2d j: %2d c6: %10.5e cn: %10.5e\n", i, j,
- C6(nbfp, atnr, i, j), C12(nbfp, atnr, i, j));
- }
- }
- }
- }
-}
-
-static void scale_box(int natoms, rvec x[], matrix box)
-{
- int i, m;
-
- if ((scale[XX] != 1.0) || (scale[YY] != 1.0) || (scale[ZZ] != 1.0))
- {
- if (debug)
- {
- fprintf(debug, "scale = %8.4f %8.4f %8.4f\n",
- scale[XX], scale[YY], scale[ZZ]);
- }
- for (m = 0; (m < DIM); m++)
- {
- box[m][m] *= scale[m];
- }
- for (i = 0; (i < natoms); i++)
- {
- for (m = 0; (m < DIM); m++)
- {
- x[i][m] *= scale[m];
- }
- }
- }
-}
-
-gmx_bool update_forcefield(FILE *fplog,
- int nfile, const t_filenm fnm[], t_forcerec *fr,
- int natoms, rvec x[], matrix box)
-{
- static int ntry, ntried;
- int i, j;
- gmx_bool bDone;
-
- /* First time around we have to read the parameters */
- if (nparm == 0)
- {
- range = read_range(ftp2fn(efDAT, nfile, fnm), &nparm);
- if (nparm == 0)
- {
- gmx_fatal(FARGS, "No correct parameter info in %s", ftp2fn(efDAT, nfile, fnm));
- }
- snew(param_val, nparm);
-
- if (opt2bSet("-ga", nfile, fnm))
- {
- /* Genetic algorithm time */
- ga = init_ga(fplog, opt2fn("-ga", nfile, fnm), nparm, range);
- }
- else
- {
- /* Determine the grid size */
- ntry = 1;
- for (i = 0; (i < nparm); i++)
- {
- ntry *= range[i].np;
- }
- ntried = 0;
-
- fprintf(fplog, "Going to try %d different combinations of %d parameters\n",
- ntry, nparm);
- }
- }
- if (ga)
- {
- update_ga(fplog, range, ga);
- }
- else
- {
- /* Increment the counter
- * Non-trivial, since this is nparm nested loops in principle
- */
- for (i = 0; (i < nparm); i++)
- {
- if (param_val[i] < (range[i].np-1))
- {
- param_val[i]++;
- for (j = 0; (j < i); j++)
- {
- param_val[j] = 0;
- }
- ntried++;
- break;
- }
- }
- if (i == nparm)
- {
- fprintf(fplog, "Finished with %d out of %d iterations\n", ntried+1, ntry);
- return TRUE;
- }
- }
-
- /* Now do the real updating */
- update_ff(fr, nparm, range, param_val);
-
- /* Update box and coordinates if necessary */
- scale_box(natoms, x, box);
-
- return FALSE;
-}
-
-static void print_range(FILE *fp, tensor P, real MSF, real energy)
-{
- int i;
-
- fprintf(fp, "%8.3f %8.3f %8.3f %8.3f",
- cost(P, MSF, energy), trace(P)/3, MSF, energy);
- for (i = 0; (i < nparm); i++)
- {
- fprintf(fp, " %s %10g", esenm[range[i].ptype], range[i].rval);
- }
- fprintf(fp, " FF\n");
- fflush(fp);
-}
-
-static real msf(int n, rvec f1[], rvec f2[])
-{
- int i, j;
- rvec ff2;
- real msf1 = 0;
-
- for (i = 0; (i < n); )
- {
- clear_rvec(ff2);
- for (j = 0; ((j < ff.molsize) && (i < n)); j++, i++)
- {
- rvec_inc(ff2, f1[i]);
- if (f2)
- {
- rvec_inc(ff2, f2[i]);
- }
- }
- msf1 += iprod(ff2, ff2);
- }
-
- return msf1/n;
-}
-
-static void print_grid(FILE *fp, real ener[], int natoms, rvec f[], rvec fshake[],
- rvec x[], t_block *mols, real mass[], tensor pres)
-{
- static gmx_bool bFirst = TRUE;
- static const char *desc[] = {
- "------------------------------------------------------------------------",
- "In the output from the forcefield scan we have the potential energy,",
- "then the root mean square force on the atoms, and finally the parameters",
- "in the order they appear in the input file.",
- "------------------------------------------------------------------------"
- };
- real msf1;
- int i;
-
- if (bFirst)
- {
- for (i = 0; (i < asize(desc)); i++)
- {
- fprintf(fp, "%s\n", desc[i]);
- }
- fflush(fp);
- bFirst = FALSE;
- }
- if ((ff.tol == 0) || (fabs(ener[F_EPOT]/ff.nmol-ff.epot) < ff.tol))
- {
- msf1 = msf(natoms, f, fshake);
- if ((ff.f_max == 0) || (msf1 < sqr(ff.f_max)))
- {
- print_range(fp, pres, msf1, ener[F_EPOT]/ff.nmol);
- }
- }
-}
-
-gmx_bool print_forcefield(FILE *fp, real ener[], int natoms, rvec f[], rvec fshake[],
- rvec x[], t_block *mols, real mass[], tensor pres)
-{
- real msf1;
-
- if (ga)
- {
- msf1 = msf(natoms, f, fshake);
- if (debug)
- {
- fprintf(fp, "Pressure: %12g, RMSF: %12g, Energy-Epot: %12g, cost: %12g\n",
- ener[F_PRES], sqrt(msf1), ener[F_EPOT]/ff.nmol-ff.epot,
- cost(pres, msf1, ener[F_EPOT]/ff.nmol));
- }
- if (print_ga(fp, ga, msf1, pres, scale, (ener[F_EPOT]/ff.nmol), range, ff.tol))
- {
- return TRUE;
- }
- fflush(fp);
- }
- else
- {
- print_grid(fp, ener, natoms, f, fshake, x, mols, mass, pres);
- }
- return FALSE;
-}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff