Port My code to Latest Gromacs

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)

Go to the main branch, commit 53deabeed3. Create your own branch from it.

Apply the patch that is attached to this message (git apply module_fhmd.dat)
module_fhmd.dat (25.9 KB)

Look at the files in src/gromacs/fhmd. This is a basic module template, upon which I added a few basic modification that you can use to port your code.

This replaces the normal integration step, by calling a function in the FHMD module. This function does nothing currently - so the simulation crashes promptly, but your code will go there. There is a bunch of printf statement starting with [EXAMPLE], which you can see when executing mdrun, and also search for in the code, in the places that might be relevant for you.

You need to toggle on the mdp parameter that enables the module for it to be active: fhmd-enabled = yes.

Replacing the update step is not something there is an official interface for, so I quickly made up an interface hack for that that works with MDModules. This means you can port your code right now without waiting for an official interface. It also mean there are known incompatibilities (see remarks below), and certainly unknown incompatibilities, and that there might occasionally be breaking changes/conflict in md.cpp that you will have to resolve, but with your code separated from md.cpp, I think it will not be too bad to keep it up to date.

Please test your port with simple cases (system with a few particle where you can indenpendently calculate what the forces and positions should be for a few step), including with/without various feature like GPU non-bonded (-pme cpu/gpu -nb cpu/gpu) that could cause trouble. In general, I give no guarantee the template code I provided works as-is, or at all.

There may be a proper interface at some point in the future - in which case the hack becomes uncessary. (You could even help develop such interface - maybe to join the weekly developper meeting sometime)

Misc remarks:

  • To see the printf when starting out, use the -ntomp 1 -ntmpi 1 command line args to get single thread execution in mdrun

  • Coordinate update on GPU is now default, but unsupported, so use also -update cpu.

  • If I’m not mistaken, the Velocity Verlet integrator is not supported with your code in the md.cpp you provided - it is also not in the template attached

  • Ultimately, automated unit test and high level test are the way to go for correctness testing.

  • Other people interested in making a module can selectively apply changes from this patch too - though they should probably not apply the changes in md.cpp and mdtype/*, and rename all references to FHMD. (Probably we should have a branch with an example module, always rebased on top of main and current release branches, for exactly this purpose of providing a quick start to new developers?)

Hi ebriand,

I do appreciate your help very much. Following your guidance, the code worked as you said. Now I am going to port my code.

I am not very familiar C++, but your example is a good starting point for me to get involved in gromacs development.

(Probably we should have a branch with an example module, always rebased on top of main and current release branches, for exactly this purpose of providing a quick start to new developers?)

In my opinion, this is a really good ideal, as it will help new developers like me very much.