Hi Gromacs developers,
I have a new particle method implemented in Gromacs 2016.4. In the new particle method, rather than updating the coordinates and velocities in accordance with conventional MD equations, the coordinates and velocities are updated following the equations (9) in this paper
https://aip.scitation.org/doi/10.1063/1.4923011
In addition, some of functions in Gromacs was invoked in my code. For example, the function to deal with PBC boundary condition is invoked many times in my code. The way used to implement the code in Gromacs 2016.4 is modified the md.cpp directly and add necessary functions to call the code. Below I showed the changes in md.cpp in Gromacs 2016.4 (in particular, the new function fhmd_update_coords to updata the coordinates and velocities) and the full md.cpp can be found in the attachment.
/*
* FHMD Initialization
*/
int is_fhmd = fhmd_init(state->box, mdatoms->homenr, mdatoms->massT, state->x, ir->delta_t, top_global, cr, &fhmd);
/* and stop now if we should */
bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
while (!bLastStep)
{
/* Determine if this is a neighbor search step */
bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
if (bPMETune && bNStList)
{
/* PME grid + cut-off optimization with GPUs or PME nodes */
pme_loadbal_do(pme_loadbal, cr,
(bVerbose && MASTER(cr)) ? stderr : NULL,
fplog,
ir, fr, state,
wcycle,
step, step_rel,
&bPMETunePrinting);
}
if(!is_fhmd)
{
update_coords(fplog, step, ir, mdatoms, state, f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
}
else /******************** FHMD coupling ********************/
{
fhmd.step_MD = step;
/*
* FHMD: Find protein COM if necessary
*/
if(fhmd.S_function == moving_sphere)
fhmd_find_protein_com(top_global, mdatoms->homenr, state->x, mdatoms->massT, cr, &fhmd);
/*
* FHMD: Estimate S in the cells and cell faces
*/
FH_S_precise(&fhmd); // FH_S(&fhmd), FH_S_weighted(&fhmd) or FH_S_precise(&fhmd)
/*
* FHMD: Update MD variables in the FH cells
*/
fhmd_update_MD_in_FH(state->x, state->v, mdatoms->massT, f, mdatoms->homenr, &fhmd);
/*
* FHMD: FH/MD Coupling
*/
if(fhmd.scheme == One_Way)
{
if(MASTER(cr) && !(fhmd.step_MD % fhmd.FH_step))
FH_do_single_timestep(&fhmd);
}
else if(fhmd.scheme == Two_Way)
{
if(PAR(cr)) fhmd_sum_arrays(cr, &fhmd);
if(fhmd.step_MD == 0)
FH_init(&fhmd, cr);
if(MASTER(cr))
{
if(fhmd.step_MD == 0)
{
FH_predictor(&fhmd);
}
else if(!(fhmd.step_MD % fhmd.FH_step))
{
FH_corrector(&fhmd);
FH_predictor(&fhmd);
}
}
}
/*
* FHMD: Broadcast new FH variables
*/
if(PAR(cr))
{
if(fhmd.scheme == One_Way)
fhmd_sum_arrays(cr, &fhmd); // for selected arrays
gmx_bcast(sizeof(FH_arrays)*fhmd.Ntot, fhmd.arr, cr); // actually we need this for the pure FH arrays only
}
/*
* FHMD: Estimate MD/FH coupling terms
*/
fhmd_calculate_MDFH_terms(&fhmd);
/*
* FHMD: Save all MD/FH arrays to files
*/
if(MASTER(cr) && (fhmd.Noutput > 0))
if(!(fhmd.step_MD % fhmd.Noutput))
{
fhmd_dump_all(&fhmd);
#ifdef FHMD_TECPLOT
fhmd_write_tecplot_data(&fhmd, fhmd.step_MD/fhmd.FH_step, (double)(fhmd.step_MD/fhmd.FH_step)*fhmd.dt_FH);
#endif
}
/*
* FHMD: Collect and print statistics
*/
fhmd_print_statistics(&fhmd, cr);
/*
* FHMD: Modified update_coords() here
*/
fhmd_update_coords(fplog, step, ir, mdatoms, state, f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr, &fhmd);
}
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
fr->bMolPBC, graph, f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
if (ir->eI == eiVVAK)
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, lastbox,
NULL, &bSumEkinhOld,
(bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again again */
copy_rvecn(cbuf, state->x, 0, state->natoms);
update_coords(fplog, step, ir, mdatoms, state, f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
The full code can be found here
https://github.com/ikorotkin/gromacs_fhmd
I am trying to port the code in latest Gromacs so that it is easy to update once the new version Gromacs is released. Meanwhile, the way to port my code should also facilitate adding a new temperature coupling method compatible with particle method. Could any one give some guidance and examples to do it?
md.log (78.1 KB)