From 685fa5f010c1419170052dff3f0a2188a77b2469 Mon Sep 17 00:00:00 2001 From: Michael Shirts Date: Fri, 1 Jan 2010 13:58:31 -0500 Subject: [PATCH] More fixes for the multicomponent lambda. --- src/gmxlib/bondfree.c | 11 ++++++++--- src/gmxlib/nonbonded/nb_free_energy.c | 8 +++++--- src/kernel/grompp.c | 28 +++++++++++++++------------ src/kernel/readir.c | 6 ++++++ src/mdlib/ns.c | 2 +- src/mdlib/sim_util.c | 36 ++++++++++++++++++++++++++--------- 6 files changed, 63 insertions(+), 28 deletions(-) diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c index 4ebe2418f9..a1d688b840 100644 --- a/src/gmxlib/bondfree.c +++ b/src/gmxlib/bondfree.c @@ -1262,10 +1262,15 @@ real posres(int nbonds, *dvdl += harmonic(pr->posres.fcA[m],pr->posres.fcB[m], 0,dpdl[m],dx[m],lambda,&v,&fm); vtot += v; - f[ai][m] += fm; - + if (f!=NULL) + { + f[ai][m] += fm; + } /* Here we correct for the pbc_dx which included rdist */ - vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm; + if (vir_diag != NULL) + { + vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm; + } } } diff --git a/src/gmxlib/nonbonded/nb_free_energy.c b/src/gmxlib/nonbonded/nb_free_energy.c index 126fead1e1..102faf6f03 100644 --- a/src/gmxlib/nonbonded/nb_free_energy.c +++ b/src/gmxlib/nonbonded/nb_free_energy.c @@ -103,7 +103,7 @@ gmx_nb_free_energy_kernel(int icoul, int n0,n1C,n1V,nnn; real Y,F,G,H,Fp,Geps,Heps2,epsC,eps2C,epsV,eps2V,VV,FF; double isp=0.564189583547756; - + bool bFEWALD; /* fix compiler warnings */ nj1 = 0; @@ -130,6 +130,8 @@ gmx_nb_free_energy_kernel(int icoul, do_tab = do_coultab || do_vdwtab; + bFEWALD = ((icoul==enbcoulFEWALD) && (alpha_coul > 0)); + /* we always use the combined table here */ tab_elemsize = 12; @@ -245,7 +247,7 @@ gmx_nb_free_energy_kernel(int icoul, n1V = tab_elemsize*n0; } - if(icoul==enbcoulOOR || icoul==enbcoulFEWALD) + if(icoul==enbcoulOOR || bFEWALD) { /* simple cutoff */ Vcoul[i] = qq[i]*rinvC; @@ -315,7 +317,7 @@ gmx_nb_free_energy_kernel(int icoul, Fscal = 0; - if (icoul==enbcoulFEWALD) + if (bFEWALD) { /* Soft-core Ewald interactions are special: * For the direct space interactions we effectively want the diff --git a/src/kernel/grompp.c b/src/kernel/grompp.c index 44e179eb5b..d6750db75a 100644 --- a/src/kernel/grompp.c +++ b/src/kernel/grompp.c @@ -1273,21 +1273,25 @@ int main (int argc, char *argv[]) } } - if (bVerbose) - fprintf(stderr,"writing run input file...\n"); - print_warn_num(TRUE); - state.fep_state = ir->fepvals->init_fep_state; - for (i=0;iefep!= efepNO) { - state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state]; - /* overwrite with explicitly stated init_lambda */ - if (ir->fepvals->init_lambda >= 0) + state.fep_state = ir->fepvals->init_fep_state; + for (i=0;ifepvals->init_lambda; - } - } - + state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state]; + /* overwrite with explicitly stated init_lambda */ + if (ir->fepvals->init_lambda >= 0) + { + state.lambda[i] = ir->fepvals->init_lambda; + } + } + } + print_warn_num(TRUE); + + if (bVerbose) + fprintf(stderr,"writing run input file...\n"); + write_tpx_state(ftp2fn(efTPX,NFILE,fnm),ir,&state,sys); thanx(stderr); diff --git a/src/kernel/readir.c b/src/kernel/readir.c index e25c5d9a90..b21228e07d 100644 --- a/src/kernel/readir.c +++ b/src/kernel/readir.c @@ -260,6 +260,11 @@ void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts, ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0)))); } } + + if (fep->bScCoul) + { + warning("Coulomb soft core not exact because of the reciprocal space calculation, as is not advised. If you wish to use, decrease the reciprocal space energy, and increase the cutoff radius."); + } /* other FEP checks that might need to be added . . . */ } @@ -1072,6 +1077,7 @@ void get_ir(const char *mdparin,const char *mdparout, RTYPE ("init-lambda", ir->fepvals->init_lambda,-1); /* start with -1 so we can recognize if it was not entered */ ITYPE ("init-fep-state", ir->fepvals->init_fep_state,0); RTYPE ("delta-lambda",ir->fepvals->delta_lambda,0.0); + /*STYPE ("foreign-lambda", fep_lambda[efptFEP], NULL);*/ /* included for backwards compatibility -- but fep-lambda overwrites */ STYPE ("fep-lambda", fep_lambda[efptFEP], NULL); STYPE ("mass-lambda", fep_lambda[efptMASS], NULL); STYPE ("coul-lambda", fep_lambda[efptCOUL], NULL); diff --git a/src/mdlib/ns.c b/src/mdlib/ns.c index 1ec716cf13..8c6d6c6d99 100644 --- a/src/mdlib/ns.c +++ b/src/mdlib/ns.c @@ -368,7 +368,7 @@ void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr) { icoulf = icoul; } - + init_nblist(&nbl->nlist_sr[eNL_VDWQQ_FREE],&nbl->nlist_lr[eNL_VDWQQ_FREE], maxsr,maxlr,ivdw,icoulf,TRUE,enlistATOM_ATOM); init_nblist(&nbl->nlist_sr[eNL_VDW_FREE],&nbl->nlist_lr[eNL_VDW_FREE], diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c index cfbbf1ddcb..d7f16b9ae1 100644 --- a/src/mdlib/sim_util.c +++ b/src/mdlib/sim_util.c @@ -406,9 +406,9 @@ void do_force(FILE *fplog,t_commrec *cr, bool bDoLongRange,bDoForces,bSepLRF; matrix boxs; real e,v,dvdl[efptNR]; + real dvdl_dum,lambda_dum; t_pbc pbc; float cycles_ppdpme,cycles_pme,cycles_seppme,cycles_force; - start = mdatoms->start; homenr = mdatoms->homenr; @@ -688,15 +688,33 @@ void do_force(FILE *fplog,t_commrec *cr, interaction_function[F_POSRES].longname,v,dvdl); } 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. - * Can we generalize this? T'would be useful. MRS. Would need to add - in the ability to read multiple position restraint files. - */ - enerd->dvdl_lin[efptRESTRAINT] += dvdl[efptRESTRAINT]; + enerd->dvdl_lin[efptRESTRAINT] += dvdl[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; in_lambda; i++) + { + if (i==0) + { + lambda_dum = lambda[efptRESTRAINT]; + } + else + { + lambda_dum = 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 */ -- 2.11.4.GIT