1patch -u -l -b -F 5 -N --suffix=.preplumed "./src/kernel/md.c" << \EOF_EOF 2--- ./src/kernel/md.c.preplumed 3+++ ./src/kernel/md.c 4@@ -86,10 +86,16 @@ 5 #include "checkpoint.h" 6 #include "mtop_util.h" 7 #include "sighandler.h" 8 #include "string2.h" 9 10+/* PLUMED */ 11+#include "../../Plumed.h" 12+extern int plumedswitch; 13+extern plumed plumedmain; 14+/* END PLUMED */ 15+ 16 #ifdef GMX_LIB_MPI 17 #include <mpi.h> 18 #endif 19 #ifdef GMX_THREADS 20 #include "tmpi.h" 21@@ -193,10 +199,16 @@ 22 #ifdef GMX_FAHCORE 23 /* Temporary addition for FAHCORE checkpointing */ 24 int chkpt_ret; 25 #endif 26 27+ /* PLUMED */ 28+ int plumedNeedsEnergy=0; 29+ int plumedWantsToStop=0; 30+ matrix plumed_vir; 31+ /* END PLUMED */ 32+ 33 /* Check for special mdrun options */ 34 bRerunMD = (Flags & MD_RERUN); 35 bIonize = (Flags & MD_IONIZE); 36 bFFscan = (Flags & MD_FFSCAN); 37 bAppend = (Flags & MD_APPENDFILES); 38@@ -440,13 +452,14 @@ 39 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy, 40 "repl_ex_nst",&repl_ex_nst); 41 /* This check needs to happen before inter-simulation 42 * signals are initialized, too */ 43 } 44- if (repl_ex_nst > 0 && MASTER(cr)) 45+ if (repl_ex_nst > 0 && MASTER(cr)){ 46 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir, 47 repl_ex_nst,repl_ex_seed); 48+ } 49 50 if (!ir->bContinuation && !bRerunMD) 51 { 52 if (mdatoms->cFREEZE && (state->flags & (1<<estV))) 53 { 54@@ -588,10 +601,60 @@ 55 } 56 } 57 fprintf(fplog,"\n"); 58 } 59 60+ /* PLUMED */ 61+ if(plumedswitch){ 62+ /* detect plumed API version */ 63+ int pversion=0; 64+ plumed_cmd(plumedmain,"getApiVersion",&pversion); 65+ /* setting kbT is only implemented with api>1) */ 66+ real kbT=ir->opts.ref_t[0]*BOLTZ; 67+ if(pversion>1) plumed_cmd(plumedmain,"setKbT",&kbT); 68+ if(pversion>2){ 69+ int res=1; 70+ if( (Flags & MD_STARTFROMCPT) ) plumed_cmd(plumedmain,"setRestart",&res); 71+ } 72+ if(cr->ms && cr->ms->nsim>1) { 73+ if(MASTER(cr)) plumed_cmd(plumedmain,"GREX setMPIIntercomm",&cr->ms->mpi_comm_masters); 74+ if(PAR(cr)){ 75+ if(DOMAINDECOMP(cr)) { 76+ plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->dd->mpi_comm_all); 77+ }else{ 78+ plumed_cmd(plumedmain,"GREX setMPIIntracomm",&cr->mpi_comm_mysim); 79+ } 80+ } 81+ plumed_cmd(plumedmain,"GREX init",NULL); 82+ } 83+ if(PAR(cr)){ 84+ if(DOMAINDECOMP(cr)) { 85+ plumed_cmd(plumedmain,"setMPIComm",&cr->dd->mpi_comm_all); 86+ }else{ 87+ plumed_cmd(plumedmain,"setMPIComm",&cr->mpi_comm_mysim); 88+ } 89+ } 90+ plumed_cmd(plumedmain,"setNatoms",&top_global->natoms); 91+ plumed_cmd(plumedmain,"setMDEngine","gromacs"); 92+ plumed_cmd(plumedmain,"setLog",fplog); 93+ real real_delta_t; 94+ real_delta_t=ir->delta_t; 95+ plumed_cmd(plumedmain,"setTimestep",&real_delta_t); 96+ plumed_cmd(plumedmain,"init",NULL); 97+ 98+ if(PAR(cr)){ 99+ if(DOMAINDECOMP(cr)) { 100+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home); 101+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex); 102+ }else{ 103+ plumed_cmd(plumedmain,"setAtomsNlocal",&mdatoms->homenr); 104+ plumed_cmd(plumedmain,"setAtomsContiguous",&mdatoms->start); 105+ } 106+ } 107+ } 108+ /* END PLUMED */ 109+ 110 /* Set and write start time */ 111 runtime_start(runtime); 112 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime); 113 wallcycle_start(wcycle,ewcRUN); 114 if (fplog) 115@@ -898,10 +961,17 @@ 116 state,&f,mdatoms,top,fr, 117 vsite,shellfc,constr, 118 nrnb,wcycle,do_verbose); 119 wallcycle_stop(wcycle,ewcDOMDEC); 120 /* If using an iterative integrator, reallocate space to match the decomposition */ 121+ 122+ /* PLUMED */ 123+ if(plumedswitch){ 124+ plumed_cmd(plumedmain,"setAtomsNlocal",&cr->dd->nat_home); 125+ plumed_cmd(plumedmain,"setAtomsGatindex",cr->dd->gatindex); 126+ } 127+ /* END PLUMED */ 128 } 129 } 130 131 if (MASTER(cr) && do_log && !bFFscan) 132 { 133@@ -1038,16 +1108,50 @@ 134 * in do_force. 135 * This is parallellized as well, and does communication too. 136 * Check comments in sim_util.c 137 */ 138 139+ 140+ /* PLUMED */ 141+ plumedNeedsEnergy=0; 142+ if(plumedswitch){ 143+ long int lstep=step; plumed_cmd(plumedmain,"setStepLong",&lstep); 144+ plumed_cmd(plumedmain,"setPositions",&state->x[mdatoms->start][0]); 145+ plumed_cmd(plumedmain,"setMasses",&mdatoms->massT[mdatoms->start]); 146+ plumed_cmd(plumedmain,"setCharges",&mdatoms->chargeA[mdatoms->start]); 147+ plumed_cmd(plumedmain,"setBox",&state->box[0][0]); 148+ plumed_cmd(plumedmain,"prepareCalc",NULL); 149+ plumed_cmd(plumedmain,"setStopFlag",&plumedWantsToStop); 150+ plumed_cmd(plumedmain,"setForces",&f[mdatoms->start][0]); 151+ plumed_cmd(plumedmain,"isEnergyNeeded",&plumedNeedsEnergy); 152+ clear_mat(plumed_vir); 153+ plumed_cmd(plumedmain,"setVirial",&plumed_vir[0][0]); 154+ } 155+ /* END PLUMED */ 156 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups, 157 state->box,state->x,&state->hist, 158 f,force_vir,mdatoms,enerd,fcd, 159 state->lambda,graph, 160 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii, 161- (bNS ? GMX_FORCE_NS : 0) | force_flags); 162+ (plumedNeedsEnergy? GMX_FORCE_ENERGY : 0) |(bNS ? GMX_FORCE_NS : 0) | force_flags); 163+ /* PLUMED */ 164+ if(plumedswitch){ 165+ if(plumedNeedsEnergy){ 166+ msmul(force_vir,2.0,plumed_vir); 167+ plumed_cmd(plumedmain,"setEnergy",&enerd->term[F_EPOT]); 168+ plumed_cmd(plumedmain,"performCalc",NULL); 169+ msmul(plumed_vir,0.5,force_vir); 170+ } else { 171+ msmul(plumed_vir,0.5,plumed_vir); 172+ plumed_cmd(plumedmain,"performCalc",NULL); 173+ m_add(force_vir,plumed_vir,force_vir); 174+ } 175+ if ((repl_ex_nst > 0) && (step > 0) && !bLastStep && 176+ do_per_step(step,repl_ex_nst)) plumed_cmd(plumedmain,"GREX savePositions",NULL); 177+ if(plumedWantsToStop) ir->nsteps=step_rel+1; 178+ } 179+ /* END PLUMED */ 180 } 181 182 GMX_BARRIER(cr->mpi_comm_mygroup); 183 184 if (bTCR) 185EOF_EOF 186patch -u -l -b -F 5 -N --suffix=.preplumed "./src/kernel/mdrun.c" << \EOF_EOF 187--- ./src/kernel/mdrun.c.preplumed 188+++ ./src/kernel/mdrun.c 189@@ -54,10 +54,16 @@ 190 #endif 191 192 /* afm stuf */ 193 #include "pull.h" 194 195+/* PLUMED */ 196+#include "../../Plumed.h" 197+int plumedswitch; 198+plumed plumedmain; 199+/* END PLUMED */ 200+ 201 int main(int argc,char *argv[]) 202 { 203 const char *desc[] = { 204 #ifdef GMX_OPENMM 205 "This is an experimental release of GROMACS for accelerated", 206@@ -376,10 +382,11 @@ 207 { efXVG, "-runav", "runaver", ffOPTWR }, 208 { efXVG, "-px", "pullx", ffOPTWR }, 209 { efXVG, "-pf", "pullf", ffOPTWR }, 210 { efMTX, "-mtx", "nm", ffOPTWR }, 211 { efNDX, "-dn", "dipole", ffOPTWR }, 212+ { efDAT, "-plumed", "plumed", ffOPTRD }, /* PLUMED */ 213 { efRND, "-multidir",NULL, ffOPTRDMULT} 214 }; 215 #define NFILE asize(fnm) 216 217 /* Command line options ! */ 218@@ -663,16 +670,52 @@ 219 220 ddxyz[XX] = (int)(realddxyz[XX] + 0.5); 221 ddxyz[YY] = (int)(realddxyz[YY] + 0.5); 222 ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5); 223 224+/* PLUMED */ 225+ plumedswitch=0; 226+ if (opt2bSet("-plumed",NFILE,fnm)) plumedswitch=1; 227+ if(plumedswitch){ 228+ int plumed_is_there=0; 229+ int real_precision=sizeof(real); 230+ real energyUnits=1.0; 231+ real lengthUnits=1.0; 232+ real timeUnits=1.0; 233+ 234+ 235+ if(!plumed_installed()){ 236+ gmx_fatal(FARGS,"Plumed is not available. Check your PLUMED_KERNEL variable."); 237+ } 238+ plumedmain=plumed_create(); 239+ plumed_cmd(plumedmain,"setRealPrecision",&real_precision); 240+// this is not necessary for gromacs units: 241+ plumed_cmd(plumedmain,"setMDEnergyUnits",&energyUnits); 242+ plumed_cmd(plumedmain,"setMDLengthUnits",&lengthUnits); 243+ plumed_cmd(plumedmain,"setMDTimeUnits",&timeUnits); 244+// 245+ plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,NFILE,fnm)); 246+ plumedswitch=1; 247+ } 248+/* END PLUMED */ 249+ 250+ 251+ 252 rc = mdrunner(nthreads, fplog,cr,NFILE,fnm,oenv,bVerbose,bCompact, 253 nstglobalcomm, ddxyz,dd_node_order,rdd,rconstr, 254 dddlb_opt[0],dlb_scale,ddcsx,ddcsy,ddcsz, 255 nstepout,resetstep,nmultisim,repl_ex_nst,repl_ex_seed, 256 pforce, cpt_period,max_hours,deviceOptions,Flags); 257 258+/* PLUMED */ 259+ if(plumedswitch){ 260+ plumed_finalize(plumedmain); 261+ } 262+/* END PLUMED */ 263+ 264+ 265+ 266 if (gmx_parallel_env_initialized()) 267 gmx_finalize(); 268 269 if (MULTIMASTER(cr)) { 270 thanx(stderr); 271EOF_EOF 272patch -u -l -b -F 5 -N --suffix=.preplumed "./src/kernel/repl_ex.c" << \EOF_EOF 273--- ./src/kernel/repl_ex.c.preplumed 274+++ ./src/kernel/repl_ex.c 275@@ -49,10 +49,16 @@ 276 #include "names.h" 277 #include "mvdata.h" 278 #include "domdec.h" 279 #include "partdec.h" 280 281+/* PLUMED */ 282+#include "../../Plumed.h" 283+extern int plumedswitch; 284+extern plumed plumedmain; 285+/* END PLUMED */ 286+ 287 typedef struct gmx_repl_ex 288 { 289 int repl; 290 int nrepl; 291 real temp; 292@@ -80,18 +86,20 @@ 293 294 snew(qall,ms->nsim); 295 qall[re->repl] = q; 296 gmx_sum_sim(ms->nsim,qall,ms); 297 298- bDiff = FALSE; 299- for(s=1; s<ms->nsim; s++) 300- { 301- if (qall[s] != qall[0]) 302- { 303+ /* PLUMED */ 304+ //for(s=1; s<ms->nsim; s++) 305+ //{ 306+ // if (qall[s] != qall[0]) 307+ // { 308 bDiff = TRUE; 309- } 310- } 311+ // } 312+ //} 313+ /* END PLUMED */ 314+ 315 if (bDiff) 316 { 317 if (re->type >= 0 && re->type < ereNR) 318 { 319 gmx_fatal(FARGS,"For replica exchange both %s and %s differ", 320@@ -230,26 +238,31 @@ 321 /* Make an index for increasing temperature order */ 322 for(i=0; i<re->nrepl; i++) 323 { 324 re->ind[i] = i; 325 } 326- for(i=0; i<re->nrepl; i++) 327- { 328- for(j=i+1; j<re->nrepl; j++) 329- { 330- if (re->q[re->ind[j]] < re->q[re->ind[i]]) 331- { 332- k = re->ind[i]; 333- re->ind[i] = re->ind[j]; 334- re->ind[j] = k; 335- } 336- else if (re->q[re->ind[j]] == re->q[re->ind[i]]) 337- { 338- gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]); 339- } 340- } 341- } 342+ /* PLUMED */ 343+ // plumed2: check if we want alternative patterns (i.e. for bias-exchange metaD) 344+ // in those cases replicas can share the same temperature. 345+ //for(i=0; i<re->nrepl; i++) 346+ // { 347+ // for(j=i+1; j<re->nrepl; j++) 348+ // { 349+ // if (re->q[re->ind[j]] < re->q[re->ind[i]]) 350+ // { 351+ // k = re->ind[i]; 352+ // re->ind[i] = re->ind[j]; 353+ // re->ind[j] = k; 354+ // } 355+ // else if (re->q[re->ind[j]] == re->q[re->ind[i]]) 356+ // { 357+ // gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]); 358+ // } 359+ // } 360+ // } 361+ //} 362+ /* END PLUMED */ 363 fprintf(fplog,"Repl "); 364 for(i=0; i<re->nrepl; i++) 365 { 366 fprintf(fplog," %3d ",re->ind[i]); 367 } 368@@ -603,10 +616,37 @@ 369 snew(bEx,re->nrepl); 370 snew(prob,re->nrepl); 371 372 exchange = -1; 373 m = (step / re->nst) % 2; 374+ /* PLUMED */ 375+ if(plumedswitch){ 376+ int partner=re->repl; 377+ int test=0; 378+ plumed_cmd(plumedmain,"getExchangesFlag",&test); 379+ if(test>0){ 380+ int *list; 381+ snew(list,re->nrepl); 382+ plumed_cmd(plumedmain,"setNumberOfReplicas",&(re->nrepl)); 383+ plumed_cmd(plumedmain,"getExchangesList",list); 384+ for(i=0; i<re->nrepl; i++) re->ind[i]=list[i]; 385+ sfree(list); 386+ } 387+ 388+ for(i=1; i<re->nrepl; i++) { 389+ if (i % 2 != m) continue; 390+ a = re->ind[i-1]; 391+ b = re->ind[i]; 392+ if(re->repl==a) partner=b; 393+ if(re->repl==b) partner=a; 394+ } 395+ plumed_cmd(plumedmain,"GREX setPartner",&partner); 396+ plumed_cmd(plumedmain,"GREX calculate",NULL); 397+ plumed_cmd(plumedmain,"GREX shareAllDeltaBias",NULL); 398+ } 399+ /* END PLUMED */ 400+ 401 for(i=1; i<re->nrepl; i++) 402 { 403 a = re->ind[i-1]; 404 b = re->ind[i]; 405 bPrint = (re->repl==a || re->repl==b); 406@@ -620,10 +660,22 @@ 407 */ 408 ediff = Epot[b] - Epot[a]; 409 betaA = 1.0/(re->q[a]*BOLTZ); 410 betaB = 1.0/(re->q[b]*BOLTZ); 411 delta = (betaA - betaB)*ediff; 412+ /* PLUMED */ 413+ if(plumedswitch){ 414+ real adb,bdb,dplumed; 415+ char buf[300]; 416+ sprintf(buf,"GREX getDeltaBias %d",a); plumed_cmd(plumedmain,buf,&adb); 417+ sprintf(buf,"GREX getDeltaBias %d",b); plumed_cmd(plumedmain,buf,&bdb); 418+ dplumed=adb*betaA+bdb*betaB; 419+ delta+=dplumed; 420+ if (bPrint) 421+ fprintf(fplog," dplumed = %10.3e d = %10.3e",dplumed,delta); 422+ } 423+ /* END PLUMED */ 424 break; 425 case ereLAMBDA: 426 /* Here we exchange based on a linear extrapolation of dV/dlambda. 427 * We would like to have the real energies 428 * from foreign lambda calculations. 429@@ -724,10 +776,14 @@ 430 gmx_multisim_t *ms; 431 int exchange=-1,shift; 432 gmx_bool bExchanged=FALSE; 433 434 ms = cr->ms; 435+ 436+ /* PLUMED */ 437+ if(plumedswitch)plumed_cmd(plumedmain,"GREX prepare",NULL); 438+ /* END PLUMED */ 439 440 if (MASTER(cr)) 441 { 442 exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box), 443 step,time); 444EOF_EOF 445