/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
*
- *
+ *
* 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.
* 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:
* GROwing Monsters And Cloning Shrimps
*/
#include "pbc.h"
#include "chargegroup.h"
#include "vec.h"
-#include "time.h"
+#include <time.h>
#include "nrnb.h"
#include "mshift.h"
#include "mdrun.h"
#include "force.h"
#include "bondf.h"
#include "pme.h"
-#include "pppm.h"
#include "disre.h"
#include "orires.h"
#include "network.h"
#include "trnio.h"
#include "xtcio.h"
#include "copyrite.h"
-
+#include "pull_rotation.h"
+#include "gmx_random.h"
#include "mpelogging.h"
#include "domdec.h"
#include "partdec.h"
#ifdef GMX_LIB_MPI
#include <mpi.h>
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "tmpi.h"
#endif
+#include "adress.h"
#include "qmmm.h"
#if 0
typedef struct gmx_timeprint {
-
+
} t_gmx_timeprint;
#endif
#ifdef HAVE_GETTIMEOFDAY
struct timeval t;
double seconds;
-
+
gettimeofday(&t,NULL);
-
+
seconds = (double) t.tv_sec + 1e-6*(double)t.tv_usec;
-
+
return seconds;
#else
double seconds;
-
+
seconds = time(NULL);
-
+
return seconds;
#endif
}
#define difftime(end,start) ((double)(end)-(double)(start))
-void print_time(FILE *out,gmx_runtime_t *runtime,gmx_large_int_t step,
+void print_time(FILE *out,gmx_runtime_t *runtime,gmx_large_int_t step,
t_inputrec *ir, t_commrec *cr)
{
time_t finish;
char timebuf[STRLEN];
double dt;
char buf[48];
-
-#ifndef GMX_THREADS
+
+#ifndef GMX_THREAD_MPI
if (!PAR(cr))
#endif
{
runtime->time_per_step = dt/(step - ir->init_step + 1);
}
dt = (ir->nsteps + ir->init_step - step)*runtime->time_per_step;
-
+
if (ir->nsteps >= 0)
{
if (dt >= 300)
- {
+ {
finish = (time_t) (runtime->last + dt);
gmx_ctime_r(&finish,timebuf,STRLEN);
sprintf(buf,"%s",timebuf);
ir->delta_t/1000*24*60*60/runtime->time_per_step);
}
}
-#ifndef GMX_THREADS
+#ifndef GMX_THREAD_MPI
if (PAR(cr))
{
fprintf(out,"\n");
fflush(out);
}
-#ifdef NO_CLOCK
+#ifdef NO_CLOCK
#define clock() -1
#endif
prev = runtime->proc;
runtime->proc = dclock();
-
+
diff = runtime->proc - prev;
#else
clock_t prev;
void runtime_end(gmx_runtime_t *runtime)
{
double now;
-
+
now = gmx_gettime();
-
+
runtime->proctime += set_proctime(runtime);
runtime->realtime = now - runtime->real;
runtime->real = now;
static void sum_forces(int start,int end,rvec f[],rvec flr[])
{
int i;
-
+
if (gmx_debug_at) {
pr_rvecs(debug,0,"fsr",f+start,end-start);
pr_rvecs(debug,0,"flr",flr+start,end-start);
rvec_inc(f[i],flr[i]);
}
-/*
+/*
* calc_f_el calculates forces due to an electric field.
*
- * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
+ * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
*
- * Et[] contains the parameters for the time dependent
- * part of the field (not yet used).
+ * Et[] contains the parameters for the time dependent
+ * part of the field (not yet used).
* Ex[] contains the parameters for
* the spatial dependent part of the field. You can have cool periodic
* fields in principle, but only a constant field is supported
- * now.
+ * now.
* The function should return the energy due to the electric field
* (if any) but for now returns 0.
*
rvec Ext;
real t0;
int i,m;
-
+
for(m=0; (m<DIM); m++)
{
if (Et[m].n > 0)
clear_mat(vir_part);
calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
-
- /* Calculate partial virial, for local atoms only, based on short range.
- * Total virial is computed in global_stat, called from do_md
+
+ /* Calculate partial virial, for local atoms only, based on short range.
+ * Total virial is computed in global_stat, called from do_md
*/
f_calc_vir(fplog,start,start+homenr,x,f,vir_part,graph,box);
inc_nrnb(nrnb,eNR_VIRIAL,homenr);
tensor vir_force,
t_mdatoms *mdatoms,
gmx_enerdata_t *enerd,t_fcdata *fcd,
- real lambda,t_graph *graph,
+ real *lambda,t_graph *graph,
t_forcerec *fr,gmx_vsite_t *vsite,rvec mu_tot,
double t,FILE *field,gmx_edsam_t ed,
gmx_bool bBornRadii,
{
int cg0,cg1,i,j;
int start,homenr;
- double mu[2*DIM];
+ double mu[2*DIM];
gmx_bool bSepDVDL,bStateChanged,bNS,bFillGrid,bCalcCGCM,bBS;
gmx_bool bDoLongRange,bDoForces,bSepLRF;
+ gmx_bool bDoAdressWF;
matrix boxs;
- real e,v,dvdl;
+ real e,v,dvdlambda[efptNR];
+ real dvdl_dum,lambda_dum;
t_pbc pbc;
float cycles_ppdpme,cycles_pme,cycles_seppme,cycles_force;
-
+
start = mdatoms->start;
homenr = mdatoms->homenr;
}
bStateChanged = (flags & GMX_FORCE_STATECHANGED);
- bNS = (flags & GMX_FORCE_NS) && (fr->bAllvsAll==FALSE);
+ bNS = (flags & GMX_FORCE_NS) && (fr->bAllvsAll==FALSE);
bFillGrid = (bNS && bStateChanged);
bCalcCGCM = (bFillGrid && !DOMAINDECOMP(cr));
bDoLongRange = (fr->bTwinRange && bNS && (flags & GMX_FORCE_DOLR));
bDoForces = (flags & GMX_FORCE_FORCES);
bSepLRF = (bDoLongRange && bDoForces && (flags & GMX_FORCE_SEPLRF));
+ /* should probably move this to the forcerec since it doesn't change */
+ bDoAdressWF = ((fr->adress_type!=eAdressOff));
if (bStateChanged)
{
update_forcerec(fplog,fr,box);
-
- /* Calculate total (local) dipole moment in a temporary common array.
+
+ /* Calculate total (local) dipole moment in a temporary common array.
* This makes it possible to sum them over nodes faster.
*/
calc_mu(start,homenr,
x,mdatoms->chargeA,mdatoms->chargeB,mdatoms->nChargePerturbed,
mu,mu+DIM);
}
-
- if (fr->ePBC != epbcNONE) {
+
+ if (fr->ePBC != epbcNONE) {
/* Compute shift vectors every step,
* because of pressure coupling or box deformation!
*/
if ((flags & GMX_FORCE_DYNAMICBOX) && bStateChanged)
calc_shifts(box,fr->shift_vec);
-
- if (bCalcCGCM) {
+
+ if (bCalcCGCM) {
put_charge_groups_in_box(fplog,cg0,cg1,fr->ePBC,box,
&(top->cgs),x,fr->cg_cm);
inc_nrnb(nrnb,eNR_CGCM,homenr);
inc_nrnb(nrnb,eNR_RESETX,cg1-cg0);
- }
+ }
else if (EI_ENERGY_MINIMIZATION(inputrec->eI) && graph) {
unshift_self(graph,box,x);
}
- }
+ }
else if (bCalcCGCM) {
calc_cgcm(fplog,cg0,cg1,&(top->cgs),x,fr->cg_cm);
inc_nrnb(nrnb,eNR_CGCM,homenr);
}
-
+
if (bCalcCGCM) {
if (PAR(cr)) {
move_cgcm(fplog,cr,fr->cg_cm);
* Since this is only implemented for domain decomposition
* and domain decomposition does not use the graph,
* we do not need to worry about shifting.
- */
+ */
wallcycle_start(wcycle,ewcPP_PMESENDX);
GMX_MPE_LOG(ev_send_coordinates_start);
}
gmx_pme_send_x(cr,bBS ? boxs : box,x,
- mdatoms->nChargePerturbed,lambda,
+ mdatoms->nChargePerturbed,lambda[efptCOUL],
( flags & GMX_FORCE_VIRIAL),step);
GMX_MPE_LOG(ev_send_coordinates_finish);
}
if (bStateChanged)
{
+
+ /* update adress weight beforehand */
+ if(bDoAdressWF)
+ {
+ /* need pbc for adress weight calculation with pbc_dx */
+ set_pbc(&pbc,inputrec->ePBC,box);
+ if(fr->adress_site == eAdressSITEcog)
+ {
+ update_adress_weights_cog(top->idef.iparams,top->idef.il,x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else if (fr->adress_site == eAdressSITEcom)
+ {
+ update_adress_weights_com(fplog,cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else if (fr->adress_site == eAdressSITEatomatom){
+ update_adress_weights_atom_per_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ else
+ {
+ update_adress_weights_atom(cg0,cg1,&(top->cgs),x,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+ }
+
for(i=0; i<2; i++)
{
for(j=0;j<DIM;j++)
for(j=0; j<DIM; j++)
{
mu_tot[j] =
- (1.0 - lambda)*fr->mu_tot[0][j] + lambda*fr->mu_tot[1][j];
+ (1.0 - lambda[efptCOUL])*fr->mu_tot[0][j] + lambda[efptCOUL]*fr->mu_tot[1][j];
}
}
if (bNS)
{
wallcycle_start(wcycle,ewcNS);
-
+
if (graph && bStateChanged)
{
/* Calculate intramolecular shift vectors to make molecules whole */
/* Do the actual neighbour searching and if twin range electrostatics
* also do the calculation of long range forces and energies.
*/
- dvdl = 0;
+ for (i=0;i<efptNR;i++) {dvdlambda[i] = 0;}
ns(fplog,fr,x,box,
groups,&(inputrec->opts),top,mdatoms,
- cr,nrnb,lambda,&dvdl,&enerd->grpp,bFillGrid,
+ cr,nrnb,lambda,dvdlambda,&enerd->grpp,bFillGrid,
bDoLongRange,bDoForces,bSepLRF ? fr->f_twin : f);
if (bSepDVDL)
{
- fprintf(fplog,sepdvdlformat,"LR non-bonded",0.0,dvdl);
+ fprintf(fplog,sepdvdlformat,"LR non-bonded",0.0,dvdlambda);
}
- enerd->dvdl_lin += dvdl;
-
+ enerd->dvdl_lin[efptVDW] += dvdlambda[efptVDW];
+ enerd->dvdl_lin[efptCOUL] += dvdlambda[efptCOUL];
+
wallcycle_stop(wcycle,ewcNS);
}
-
- if (inputrec->implicit_solvent && bNS)
+
+ if (inputrec->implicit_solvent && bNS)
{
make_gb_nblist(cr,inputrec->gb_algorithm,inputrec->rlist,
x,box,fr,&top->idef,graph,fr->born);
}
-
+
if (DOMAINDECOMP(cr))
{
if (!(cr->duty & DUTY_PME))
dd_force_flop_start(cr->dd,nrnb);
}
}
-
+
+ if (inputrec->bRot)
+ {
+ /* Enforced rotation has its own cycle counter that starts after the collective
+ * coordinates have been communicated. It is added to ddCyclF to allow
+ * for proper load-balancing */
+ wallcycle_start(wcycle,ewcROT);
+ do_rotation(cr,inputrec,box,x,t,step,wcycle,bNS);
+ wallcycle_stop(wcycle,ewcROT);
+ }
+
/* Start the force cycle counter.
* This counter is stopped in do_forcelow_level.
* No parallel communication should occur while this counter is running,
* since that will interfere with the dynamic load balancing.
*/
wallcycle_start(wcycle,ewcFORCE);
-
+
if (bDoForces)
{
/* Reset forces for which the virial is calculated separately:
* PME/Ewald forces if necessary */
- if (fr->bF_NoVirSum)
+ if (fr->bF_NoVirSum)
{
if (flags & GMX_FORCE_VIRIAL)
{
if ((flags & GMX_FORCE_BONDED) && top->idef.il[F_POSRES].nr > 0)
{
- /* Position restraints always require full pbc */
- set_pbc(&pbc,inputrec->ePBC,box);
+ /* Position restraints always require full pbc. Check if we already did it for Adress */
+ if(!(bStateChanged && bDoAdressWF))
+ {
+ set_pbc(&pbc,inputrec->ePBC,box);
+ }
v = posres(top->idef.il[F_POSRES].nr,top->idef.il[F_POSRES].iatoms,
top->idef.iparams_posres,
(const rvec*)x,fr->f_novirsum,fr->vir_diag_posres,
- inputrec->ePBC==epbcNONE ? NULL : &pbc,lambda,&dvdl,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc,lambda[efptRESTRAINT],&(dvdlambda[efptRESTRAINT]),
fr->rc_scaling,fr->ePBC,fr->posres_com,fr->posres_comB);
if (bSepDVDL)
{
fprintf(fplog,sepdvdlformat,
- interaction_function[F_POSRES].longname,v,dvdl);
+ interaction_function[F_POSRES].longname,v,dvdlambda);
}
enerd->term[F_POSRES] += v;
/* This linear lambda dependence assumption is only correct
* when only k depends on lambda,
* not when the reference position depends on lambda.
- * grompp checks for this.
+ * grompp checks for this. (verify this is still the case?)
*/
- enerd->dvdl_lin += dvdl;
+ enerd->dvdl_nonlin[efptRESTRAINT] += dvdlambda[efptRESTRAINT]; /* if just the force constant changes, this is linear,
+ but we can't be sure w/o additional checking that is
+ hard to do at this level of code. Otherwise,
+ the dvdl is not differentiable */
inc_nrnb(nrnb,eNR_POSRES,top->idef.il[F_POSRES].nr/2);
- }
+ if ((inputrec->fepvals->n_lambda > 0) && (flags & GMX_FORCE_DHDL))
+ {
+ for(i=0; i<enerd->n_lambda; i++)
+ {
+ lambda_dum = (i==0 ? lambda[efptRESTRAINT] : inputrec->fepvals->all_lambda[efptRESTRAINT][i-1]);
+ v = posres(top->idef.il[F_POSRES].nr,top->idef.il[F_POSRES].iatoms,
+ top->idef.iparams_posres,
+ (const rvec*)x,NULL,NULL,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc,lambda_dum,&dvdl_dum,
+ fr->rc_scaling,fr->ePBC,fr->posres_com,fr->posres_comB);
+ enerd->enerpart_lambda[i] += v;
+ }
+ }
+ }
- /* Compute the bonded and non-bonded energies and optionally forces */
+ /* Compute the bonded and non-bonded energies and optionally forces */
do_force_lowlevel(fplog,step,fr,inputrec,&(top->idef),
cr,nrnb,wcycle,mdatoms,&(inputrec->opts),
x,hist,f,enerd,fcd,mtop,top,fr->born,
&(top->atomtypes),bBornRadii,box,
- lambda,graph,&(top->excls),fr->mu_tot,
+ inputrec->fepvals,lambda,graph,&(top->excls),fr->mu_tot,
flags,&cycles_pme);
-
+
cycles_force = wallcycle_stop(wcycle,ewcFORCE);
GMX_BARRIER(cr->mpi_comm_mygroup);
-
+
if (ed)
{
do_flood(fplog,cr,x,f,ed,box,step,bNS);
}
-
+
if (DOMAINDECOMP(cr))
{
dd_force_flop_stop(cr->dd,nrnb);
dd_cycles_add(cr->dd,cycles_force-cycles_pme,ddCyclF);
}
}
-
+
if (bDoForces)
{
if (IR_ELEC_FIELD(*inputrec))
start,homenr,mdatoms->chargeA,x,fr->f_novirsum,
inputrec->ex,inputrec->et,t);
}
-
+
+ if (bDoAdressWF && fr->adress_icor == eAdressICThermoForce)
+ {
+ /* Compute thermodynamic force in hybrid AdResS region */
+ adress_thermo_force(start,homenr,&(top->cgs),x,fr->f_novirsum,fr,mdatoms,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc);
+ }
+
/* Communicate the forces */
if (PAR(cr))
{
wallcycle_stop(wcycle,ewcVSITESPREAD);
}
}
-
+
if (flags & GMX_FORCE_VIRIAL)
{
/* Calculation of the virial must be done after vsites! */
}
}
+ enerd->term[F_COM_PULL] = 0;
if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
{
/* Calculate the center of mass forces, this requires communication,
* which is why we call pull_potential after calc_virial.
*/
set_pbc(&pbc,inputrec->ePBC,box);
- dvdl = 0;
- enerd->term[F_COM_PULL] =
+ dvdlambda[efptRESTRAINT] = 0;
+ enerd->term[F_COM_PULL] +=
pull_potential(inputrec->ePull,inputrec->pull,mdatoms,&pbc,
- cr,t,lambda,x,f,vir_force,&dvdl);
+ cr,t,lambda[efptRESTRAINT],x,f,vir_force,&(dvdlambda[efptRESTRAINT]));
if (bSepDVDL)
{
- fprintf(fplog,sepdvdlformat,"Com pull",enerd->term[F_COM_PULL],dvdl);
+ fprintf(fplog,sepdvdlformat,"Com pull",enerd->term[F_COM_PULL],dvdlambda[efptRESTRAINT]);
}
- enerd->dvdl_lin += dvdl;
+ enerd->dvdl_lin[efptRESTRAINT] += dvdlambda[efptRESTRAINT];
+ }
+
+ /* Add the forces from enforced rotation potentials (if any) */
+ if (inputrec->bRot)
+ {
+ wallcycle_start(wcycle,ewcROTadd);
+ enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr,step,t);
+ wallcycle_stop(wcycle,ewcROTadd);
}
if (PAR(cr) && !(cr->duty & DUTY_PME))
cycles_ppdpme = wallcycle_stop(wcycle,ewcPPDURINGPME);
dd_cycles_add(cr->dd,cycles_ppdpme,ddCyclPPduringPME);
- /* In case of node-splitting, the PP nodes receive the long-range
+ /* In case of node-splitting, the PP nodes receive the long-range
* forces, virial and energy from the PME nodes here.
- */
+ */
wallcycle_start(wcycle,ewcPP_PMEWAITRECVF);
- dvdl = 0;
- gmx_pme_receive_f(cr,fr->f_novirsum,fr->vir_el_recip,&e,&dvdl,
+ dvdlambda[efptCOUL] = 0;
+ gmx_pme_receive_f(cr,fr->f_novirsum,fr->vir_el_recip,&e,&dvdlambda[efptCOUL],
&cycles_seppme);
if (bSepDVDL)
{
- fprintf(fplog,sepdvdlformat,"PME mesh",e,dvdl);
+ fprintf(fplog,sepdvdlformat,"PME mesh",e,dvdlambda[efptCOUL]);
}
enerd->term[F_COUL_RECIP] += e;
- enerd->dvdl_lin += dvdl;
+ enerd->dvdl_lin[efptCOUL] += dvdlambda[efptCOUL];
if (wcycle)
{
dd_cycles_add(cr->dd,cycles_seppme,ddCyclPME);
{
if (vsite)
{
- /* Spread the mesh force on virtual sites to the other particles...
+ /* Spread the mesh force on virtual sites to the other particles...
* This is parallellized. MPI communication is performed
* if the constructing atoms aren't local.
*/
}
}
}
-
+
/* Sum the potential energy terms from group contributions */
sum_epot(&(inputrec->opts),enerd);
-
+
if (fr->print_force >= 0 && bDoForces)
{
print_large_forces(stderr,mdatoms,cr,step,fr->print_force,x,f);
int i,m,start,end;
gmx_large_int_t step;
real dt=ir->delta_t;
- real dvdlambda;
+ real dvdl_dum;
rvec *savex;
-
+
snew(savex,state->natoms);
start = md->start;
end = md->homenr + start;
-
+
if (debug)
fprintf(debug,"vcm: start=%d, homenr=%d, end=%d\n",
start,md->homenr,end);
fprintf(fplog,"\nConstraining the starting coordinates (step %s)\n",
gmx_step_str(step,buf));
}
- dvdlambda = 0;
-
+ dvdl_dum = 0;
+
/* constrain the current position */
constrain(NULL,TRUE,FALSE,constr,&(top->idef),
ir,NULL,cr,step,0,md,
state->x,state->x,NULL,
- state->box,state->lambda,&dvdlambda,
+ state->box,state->lambda[efptBONDED],&dvdl_dum,
NULL,NULL,nrnb,econqCoord,ir->epc==epcMTTK,state->veta,state->veta);
- if (EI_VV(ir->eI))
+ if (EI_VV(ir->eI))
{
/* constrain the inital velocity, and save it */
/* also may be useful if we need the ekin from the halfstep for velocity verlet */
constrain(NULL,TRUE,FALSE,constr,&(top->idef),
ir,NULL,cr,step,0,md,
state->x,state->v,state->v,
- state->box,state->lambda,&dvdlambda,
+ state->box,state->lambda[efptBONDED],&dvdl_dum,
NULL,NULL,nrnb,econqVeloc,ir->epc==epcMTTK,state->veta,state->veta);
}
/* constrain the inital velocities at t-dt/2 */
if (EI_STATE_VELOCITY(ir->eI) && ir->eI!=eiVV)
{
- for(i=start; (i<end); i++)
+ for(i=start; (i<end); i++)
{
- for(m=0; (m<DIM); m++)
+ for(m=0; (m<DIM); m++)
{
/* Reverse the velocity */
state->v[i][m] = -state->v[i][m];
savex[i][m] = state->x[i][m] + dt*state->v[i][m];
}
}
- /* Shake the positions at t=-dt with the positions at t=0
- * as reference coordinates.
+ /* Shake the positions at t=-dt with the positions at t=0
+ * as reference coordinates.
*/
if (fplog)
{
fprintf(fplog,"\nConstraining the coordinates at t0-dt (step %s)\n",
gmx_step_str(step,buf));
}
- dvdlambda = 0;
+ dvdl_dum = 0;
constrain(NULL,TRUE,FALSE,constr,&(top->idef),
ir,NULL,cr,step,-1,md,
state->x,savex,NULL,
- state->box,state->lambda,&dvdlambda,
+ state->box,state->lambda[efptBONDED],&dvdl_dum,
state->v,NULL,nrnb,econqCoord,ir->epc==epcMTTK,state->veta,state->veta);
-
+
for(i=start; i<end; i++) {
for(m=0; m<DIM; m++) {
/* Re-reverse the velocities */
}
}
}
-
sfree(savex);
}
double r0,r1,r,rc3,rc9,ea,eb,ec,pa,pb,pc,pd;
double invscale,invscale2,invscale3;
int ri0,ri1,ri,i,offstart,offset;
- real scale,*vdwtab;
+ real scale,*vdwtab;
fr->enershiftsix = 0;
fr->enershifttwelve = 0;
*/
eners[0] += 4.0*M_PI*fr->enershiftsix*rc3/3.0;
eners[1] += 4.0*M_PI*fr->enershifttwelve*rc3/3.0;
-
- invscale = 1.0/(scale);
+
+ invscale = 1.0/(scale);
invscale2 = invscale*invscale;
invscale3 = invscale*invscale2;
switched function. We perform both the pressure and energy
loops at the same time for simplicity, as the computational
cost is low. */
-
+
for (i=0;i<2;i++) {
enersum = 0.0; virsum = 0.0;
if (i==0)
ea = invscale3;
eb = 2.0*invscale2*r;
ec = invscale*r*r;
-
+
pa = invscale3;
pb = 3.0*invscale2*r;
pc = 3.0*invscale*r*r;
pd = r*r*r;
-
+
/* this "8" is from the packing in the vdwtab array - perhaps
should be #define'ed? */
offset = 8*ri + offstart;
f = vdwtab[offset+1];
g = vdwtab[offset+2];
h = vdwtab[offset+3];
-
+
enersum += y0*(ea/3 + eb/2 + ec) + f*(ea/4 + eb/3 + ec/2)+
- g*(ea/5 + eb/4 + ec/3) + h*(ea/6 + eb/5 + ec/4);
- virsum += f*(pa/4 + pb/3 + pc/2 + pd) +
+ g*(ea/5 + eb/4 + ec/3) + h*(ea/6 + eb/5 + ec/4);
+ virsum += f*(pa/4 + pb/3 + pc/2 + pd) +
2*g*(pa/5 + pb/4 + pc/3 + pd/2) + 3*h*(pa/6 + pb/5 + pc/4 + pd/3);
-
+
}
enersum *= 4.0*M_PI;
- virsum *= 4.0*M_PI;
+ virsum *= 4.0*M_PI;
eners[i] -= enersum;
virs[i] -= virsum;
}
eners[1] += 4.0*M_PI/(9.0*rc9);
virs[0] += 8.0*M_PI/rc3;
virs[1] += -16.0*M_PI/(3.0*rc9);
- }
+ }
else if ((fr->vdwtype == evdwCUT) || (fr->vdwtype == evdwUSER)) {
if (fr->vdwtype == evdwUSER && fplog)
fprintf(fplog,
gmx_bool bCorrAll,bCorrPres;
real dvdlambda,invvol,dens,ninter,avcsix,avctwelve,enerdiff,svir=0,spres=0;
int m;
-
+
*prescorr = 0;
*enercorr = 0;
*dvdlcorr = 0;
-
+
clear_mat(virial);
clear_mat(pres);
-
+
if (ir->eDispCorr != edispcNO) {
bCorrAll = (ir->eDispCorr == edispcAllEner ||
ir->eDispCorr == edispcAllEnerPres);
bCorrPres = (ir->eDispCorr == edispcEnerPres ||
ir->eDispCorr == edispcAllEnerPres);
-
+
invvol = 1/det(box);
- if (fr->n_tpi)
+ if (fr->n_tpi)
{
/* Only correct for the interactions with the inserted molecule */
dens = (natoms - fr->n_tpi)*invvol;
ninter = fr->n_tpi;
- }
- else
+ }
+ else
{
dens = natoms*invvol;
ninter = 0.5*natoms;
}
-
- if (ir->efep == efepNO)
+
+ if (ir->efep == efepNO)
{
avcsix = fr->avcsix[0];
avctwelve = fr->avctwelve[0];
- }
- else
+ }
+ else
{
avcsix = (1 - lambda)*fr->avcsix[0] + lambda*fr->avcsix[1];
avctwelve = (1 - lambda)*fr->avctwelve[0] + lambda*fr->avctwelve[1];
}
-
+
enerdiff = ninter*(dens*fr->enerdiffsix - fr->enershiftsix);
*enercorr += avcsix*enerdiff;
dvdlambda = 0.0;
- if (ir->efep != efepNO)
+ if (ir->efep != efepNO)
{
dvdlambda += (fr->avcsix[1] - fr->avcsix[0])*enerdiff;
}
- if (bCorrAll)
+ if (bCorrAll)
{
enerdiff = ninter*(dens*fr->enerdifftwelve - fr->enershifttwelve);
*enercorr += avctwelve*enerdiff;
- if (fr->efep != efepNO)
+ if (fr->efep != efepNO)
{
dvdlambda += (fr->avctwelve[1] - fr->avctwelve[0])*enerdiff;
}
}
-
- if (bCorrPres)
+
+ if (bCorrPres)
{
svir = ninter*dens*avcsix*fr->virdiffsix/3.0;
if (ir->eDispCorr == edispcAllEnerPres)
}
/* The factor 2 is because of the Gromacs virial definition */
spres = -2.0*invvol*svir*PRESFAC;
-
+
for(m=0; m<DIM; m++) {
virial[m][m] += svir;
pres[m][m] += spres;
}
*prescorr += spres;
}
-
+
/* Can't currently control when it prints, for now, just print when degugging */
if (debug)
{
fprintf(debug,"Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
avcsix,avctwelve);
}
- if (bCorrPres)
+ if (bCorrPres)
{
fprintf(debug,
"Long Range LJ corr.: Epot %10g, Pres: %10g, Vir: %10g\n",
fprintf(debug,"Long Range LJ corr.: Epot %10g\n",*enercorr);
}
}
-
+
if (fr->bSepDVDL && do_per_step(step,ir->nstlog))
{
fprintf(fplog,sepdvdlformat,"Dispersion correction",
*enercorr,dvdlambda);
}
- if (fr->efep != efepNO)
+ if (fr->efep != efepNO)
{
*dvdlcorr += dvdlambda;
}
as = 0;
for(mb=0; mb<mtop->nmolblock; mb++) {
molb = &mtop->molblock[mb];
- if (molb->natoms_mol == 1 ||
+ if (molb->natoms_mol == 1 ||
(!bFirst && mtop->moltype[molb->type].cgs.nr == 1)) {
/* Just one atom or charge group in the molecule, no PBC required */
as += molb->nmol*molb->natoms_mol;
/* Pass NULL iso fplog to avoid graph prints for each molecule type */
mk_graph_ilist(NULL,mtop->moltype[molb->type].ilist,
0,molb->natoms_mol,FALSE,FALSE,graph);
-
+
for(mol=0; mol<molb->nmol; mol++) {
mk_mshift(fplog,graph,ePBC,box,x+as);
-
+
shift_self(graph,box,x+as);
/* The molecule is whole now.
* We don't need the second mk_mshift call as in do_pbc_first,
* since we no longer need this graph.
*/
-
+
as += molb->natoms_mol;
}
done_graph(graph);
#ifdef GMX_MPI
MPI_Reduce(nrnb->n,nrnb_tot->n,eNRNB,MPI_DOUBLE,MPI_SUM,
MASTERRANK(cr),cr->mpi_comm_mysim);
-#endif
+#endif
} else {
nrnb_tot = nrnb;
}
-
+
if (SIMMASTER(cr)) {
print_flop(fplog,nrnb_tot,&nbfs,&mflop);
if (cr->nnodes > 1) {
cr->mpi_comm_mysim);
}
}
-#endif
+#endif
if (SIMMASTER(cr)) {
wallcycle_print(fplog,cr->nnodes,cr->npmenodes,runtime->realtime,
} else {
delta_t = 0;
}
-
+
if (fplog) {
print_perf(fplog,runtime->proctime,runtime->realtime,
cr->nnodes-cr->npmenodes,
}
}
+extern void initialize_lambdas(FILE *fplog,t_inputrec *ir,int *fep_state,real *lambda,double *lam0)
+{
+ /* this function works, but could probably use a logic rewrite to keep all the different
+ types of efep straight. */
+
+ int i;
+ t_lambda *fep = ir->fepvals;
+
+ if ((ir->efep==efepNO) && (ir->bSimTemp == FALSE)) {
+ for (i=0;i<efptNR;i++) {
+ lambda[i] = 0.0;
+ if (lam0)
+ {
+ lam0[i] = 0.0;
+ }
+ }
+ return;
+ } else {
+ *fep_state = fep->init_fep_state; /* this might overwrite the checkpoint
+ if checkpoint is set -- a kludge is in for now
+ to prevent this.*/
+ for (i=0;i<efptNR;i++)
+ {
+ /* overwrite lambda state with init_lambda for now for backwards compatibility */
+ if (fep->init_lambda>=0) /* if it's -1, it was never initializd */
+ {
+ lambda[i] = fep->init_lambda;
+ if (lam0) {
+ lam0[i] = lambda[i];
+ }
+ }
+ else
+ {
+ lambda[i] = fep->all_lambda[i][*fep_state];
+ if (lam0) {
+ lam0[i] = lambda[i];
+ }
+ }
+ }
+ if (ir->bSimTemp) {
+ /* need to rescale control temperatures to match current state */
+ for (i=0;i<ir->opts.ngtc;i++) {
+ if (ir->opts.ref_t[i] > 0) {
+ ir->opts.ref_t[i] = ir->simtempvals->temperatures[*fep_state];
+ }
+ }
+ }
+ }
+
+ /* Send to the log the information on the current lambdas */
+ if (fplog != NULL)
+ {
+ fprintf(fplog,"Initial vector of lambda components:[ ");
+ for (i=0;i<efptNR;i++)
+ {
+ fprintf(fplog,"%10.4f ",lambda[i]);
+ }
+ fprintf(fplog,"]\n");
+ }
+ return;
+}
+
+
void init_md(FILE *fplog,
t_commrec *cr,t_inputrec *ir,const output_env_t oenv,
double *t,double *t0,
- real *lambda,double *lam0,
+ real *lambda, int *fep_state, double *lam0,
t_nrnb *nrnb,gmx_mtop_t *mtop,
gmx_update_t *upd,
int nfile,const t_filenm fnm[],
{
int i,j,n;
real tmpt,mod;
-
+
/* Initial values */
*t = *t0 = ir->init_t;
- if (ir->efep != efepNO)
- {
- *lam0 = ir->init_lambda;
- *lambda = *lam0 + ir->init_step*ir->delta_lambda;
- }
- else
- {
- *lambda = *lam0 = 0.0;
- }
*bSimAnn=FALSE;
for(i=0;i<ir->opts.ngtc;i++)
{
update_annealing_target_temp(&(ir->opts),ir->init_t);
}
-
+
+ /* Initialize lambda variables */
+ initialize_lambdas(fplog,ir,fep_state,lambda,lam0);
+
if (upd)
{
*upd = init_update(fplog,ir);
}
-
+
+
if (vcm != NULL)
{
*vcm = init_vcm(fplog,&mtop->groups,ir);
}
-
+
if (EI_DYNAMICS(ir->eI) && !(Flags & MD_APPENDFILES))
{
if (ir->etc == etcBERENDSEN)
please_cite(fplog,"Bussi2007a");
}
}
-
+
init_nrnb(nrnb);
-
+
if (nfile != -1)
{
*outf = init_mdoutf(nfile,fnm,Flags,cr,ir,oenv);
*mdebin = init_mdebin((Flags & MD_APPENDFILES) ? NULL : (*outf)->fp_ene,
mtop,ir, (*outf)->fp_dhdl);
}
-
- /* Initiate variables */
+
+ if (ir->bAdress)
+ {
+ please_cite(fplog,"Fritsch12");
+ please_cite(fplog,"Junghans10");
+ }
+ /* Initiate variables */
clear_mat(force_vir);
clear_mat(shake_vir);
clear_rvec(mu_tot);
-
+
debug_gmx();
}
+
+
+
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff