From: Berk Hess Date: Thu, 12 Jun 2014 08:23:27 +0000 (+0200) Subject: Made shells work with the Verlet scheme X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=debd482de1251382507630e42c74898c1dff53da;p=alexxy%2Fgromacs.git Made shells work with the Verlet scheme The issue was that atoms iso charge group should be put in the box. Additionally water_pol only worked within a charge group, now fixed. Note that with shells working across charge groups, which is always the case with the Verlet scheme, no shell prediction is done. We should implement prediction, since it improves performance. Fixes #1429. Change-Id: I2ebfc2d91fc161167f8f2573b61e1f519cf11fd8 --- diff --git a/include/shellfc.h b/include/shellfc.h index 8c11aedb15..87e291d2fa 100644 --- a/include/shellfc.h +++ b/include/shellfc.h @@ -48,7 +48,6 @@ extern "C" { */ GMX_LIBMD_EXPORT gmx_shellfc_t init_shell_flexcon(FILE *fplog, - gmx_bool bCutoffSchemeIsVerlet, gmx_mtop_t *mtop, int nflexcon, rvec *x); diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c index 66a91d4ac0..2c25351ec2 100644 --- a/src/gmxlib/bondfree.c +++ b/src/gmxlib/bondfree.c @@ -736,7 +736,8 @@ real water_pol(int nbonds, * a shell connected to a dummy with spring constant that differ in the * three spatial dimensions in the molecular frame. */ - int i, m, aO, aH1, aH2, aD, aS, type, type0; + int i, m, aO, aH1, aH2, aD, aS, type, type0, ki; + ivec dt; rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj; #ifdef DEBUG rvec df; @@ -779,11 +780,11 @@ real water_pol(int nbonds, aS = forceatoms[i+5]; /* Compute vectors describing the water frame */ - rvec_sub(x[aH1], x[aO], dOH1); - rvec_sub(x[aH2], x[aO], dOH2); - rvec_sub(x[aH2], x[aH1], dHH); - rvec_sub(x[aD], x[aO], dOD); - rvec_sub(x[aS], x[aD], dDS); + pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1); + pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2); + pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH); + pbc_rvec_sub(pbc, x[aD], x[aO], dOD); + ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS); cprod(dOH1, dOH2, nW); /* Compute inverse length of normal vector @@ -838,6 +839,13 @@ real water_pol(int nbonds, kdx[YY] = kk[YY]*dx[YY]; kdx[ZZ] = kk[ZZ]*dx[ZZ]; vtot += iprod(dx, kdx); + + if (g) + { + ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt); + ki = IVEC2IS(dt); + } + for (m = 0; (m < DIM); m++) { /* This is a tensor operation but written out for speed */ @@ -848,8 +856,10 @@ real water_pol(int nbonds, #ifdef DEBUG df[m] = fij; #endif - f[aS][m] += fij; - f[aD][m] -= fij; + f[aS][m] += fij; + f[aD][m] -= fij; + fshift[ki][m] += fij; + fshift[CENTRAL][m] -= fij; } #ifdef DEBUG if (debug) diff --git a/src/kernel/md.c b/src/kernel/md.c index da68da10f2..8909fb8ed7 100644 --- a/src/kernel/md.c +++ b/src/kernel/md.c @@ -342,7 +342,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], debug_gmx(); /* Check for polarizable models and flexible constraints */ - shellfc = init_shell_flexcon(fplog, fr->cutoff_scheme == ecutsVERLET, + shellfc = init_shell_flexcon(fplog, top_global, n_flexible_constraints(constr), (ir->bContinuation || (DOMAINDECOMP(cr) && !MASTER(cr))) ? diff --git a/src/mdlib/shellfc.c b/src/mdlib/shellfc.c index 68e6c3be5e..0acd975302 100644 --- a/src/mdlib/shellfc.c +++ b/src/mdlib/shellfc.c @@ -231,7 +231,6 @@ static void predict_shells(FILE *fplog, rvec x[], rvec v[], real dt, } gmx_shellfc_t init_shell_flexcon(FILE *fplog, - gmx_bool bCutoffSchemeIsVerlet, gmx_mtop_t *mtop, int nflexcon, rvec *x) { @@ -285,11 +284,6 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog, return NULL; } - if (bCutoffSchemeIsVerlet) - { - gmx_fatal(FARGS, "The shell code does not work with the Verlet cut-off scheme.\n"); - } - snew(shfc, 1); shfc->nflexcon = nflexcon; @@ -533,6 +527,12 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog, { fprintf(fplog, "\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n"); } + /* Prediction improves performance, so we should implement either: + * 1. communication for the atoms needed for prediction + * 2. prediction using the velocities of shells; currently the + * shell velocities are zeroed, it's a bit tricky to keep + * track of the shell displacements and thus the velocity. + */ shfc->bPredict = FALSE; } } @@ -1001,34 +1001,31 @@ int relax_shell_flexcon(FILE *fplog, t_commrec *cr, gmx_bool bVerbose, force[i] = shfc->f[i]; } - /* With particle decomposition this code only works - * when all particles involved with each shell are in the same cg. - */ - if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr)) { /* This is the only time where the coordinates are used * before do_force is called, which normally puts all * charge groups in the box. */ - if (PARTDECOMP(cr)) + if (inputrec->cutoff_scheme == ecutsVERLET) { - pd_cg_range(cr, &cg0, &cg1); + put_atoms_in_box_omp(fr->ePBC, state->box, md->homenr, state->x); } else { cg0 = 0; cg1 = top->cgs.nr; + put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box, + &(top->cgs), state->x, fr->cg_cm); } - put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box, - &(top->cgs), state->x, fr->cg_cm); + if (graph) { mk_mshift(fplog, graph, fr->ePBC, state->box, state->x); } } - /* After this all coordinate arrays will contain whole molecules */ + /* After this all coordinate arrays will contain whole charge groups */ if (graph) { shift_self(graph, state->box, state->x);