#include "vsite.h"
#include "force.h"
#include "mdrun.h"
+#include "md_support.h"
#include "domdec.h"
#include "partdec.h"
#include "trnio.h"
#include "mtop_util.h"
#include "gmxfio.h"
#include "pme.h"
+#include "bondf.h"
+#include "gmx_omp_nthreads.h"
+
typedef struct {
t_state s;
}
*f_global = *f;
- if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
+ forcerec_set_excl_load(fr,*top,cr);
+
+ init_bonded_thread_force_reduction(fr,&(*top)->idef);
+
+ if (ir->ePBC != epbcNONE && !fr->bMolPBC)
{
*graph = mk_graph(fplog,&((*top)->idef),0,top_global->natoms,FALSE,FALSE);
}
dvdlambda=0;
constrain(PAR(cr) ? NULL : fplog,TRUE,TRUE,constr,&(*top)->idef,
ir,NULL,cr,-1,0,mdatoms,
- ems->s.x,ems->s.x,NULL,ems->s.box,
+ ems->s.x,ems->s.x,NULL,fr->bMolPBC,ems->s.box,
ems->s.lambda[efptFEP],&dvdlambda,
NULL,NULL,nrnb,econqCoord,FALSE,0,0);
}
{
if (!(cr->duty & DUTY_PME)) {
/* Tell the PME only node to finish */
- gmx_pme_finish(cr);
+ gmx_pme_send_finish(cr);
}
done_mdoutf(outf);
}
static void do_em_step(t_commrec *cr,t_inputrec *ir,t_mdatoms *md,
- em_state_t *ems1,real a,rvec *f,em_state_t *ems2,
- gmx_constr_t constr,gmx_localtop_t *top,
- t_nrnb *nrnb,gmx_wallcycle_t wcycle,
- gmx_large_int_t count)
+ gmx_bool bMolPBC,
+ em_state_t *ems1,real a,rvec *f,em_state_t *ems2,
+ gmx_constr_t constr,gmx_localtop_t *top,
+ t_nrnb *nrnb,gmx_wallcycle_t wcycle,
+ gmx_large_int_t count)
{
- t_state *s1,*s2;
- int start,end,gf,i,m;
- rvec *x1,*x2;
- real dvdlambda;
-
- s1 = &ems1->s;
- s2 = &ems2->s;
+ t_state *s1,*s2;
+ int i;
+ int start,end;
+ rvec *x1,*x2;
+ real dvdlambda;
- if (DOMAINDECOMP(cr) && s1->ddp_count != cr->dd->ddp_count)
- gmx_incons("state mismatch in do_em_step");
+ s1 = &ems1->s;
+ s2 = &ems2->s;
- s2->flags = s1->flags;
+ if (DOMAINDECOMP(cr) && s1->ddp_count != cr->dd->ddp_count)
+ {
+ gmx_incons("state mismatch in do_em_step");
+ }
- if (s2->nalloc != s1->nalloc) {
- s2->nalloc = s1->nalloc;
- srenew(s2->x,s1->nalloc);
- srenew(ems2->f, s1->nalloc);
- if (s2->flags & (1<<estCGP))
- srenew(s2->cg_p, s1->nalloc);
- }
+ s2->flags = s1->flags;
- s2->natoms = s1->natoms;
- /* Copy free energy state -> is this necessary? */
- for (i=0;i<efptNR;i++)
- {
- s2->lambda[i] = s1->lambda[i];
- }
- copy_mat(s1->box,s2->box);
+ if (s2->nalloc != s1->nalloc)
+ {
+ s2->nalloc = s1->nalloc;
+ srenew(s2->x,s1->nalloc);
+ srenew(ems2->f, s1->nalloc);
+ if (s2->flags & (1<<estCGP))
+ {
+ srenew(s2->cg_p, s1->nalloc);
+ }
+ }
+
+ s2->natoms = s1->natoms;
+ copy_mat(s1->box,s2->box);
+ /* Copy free energy state */
+ for (i=0;i<efptNR;i++)
+ {
+ s2->lambda[i] = s1->lambda[i];
+ }
+ copy_mat(s1->box,s2->box);
- start = md->start;
- end = md->start + md->homenr;
+ start = md->start;
+ end = md->start + md->homenr;
- x1 = s1->x;
- x2 = s2->x;
- gf = 0;
- for(i=start; i<end; i++) {
- if (md->cFREEZE)
- gf = md->cFREEZE[i];
- for(m=0; m<DIM; m++) {
- if (ir->opts.nFreeze[gf][m])
- x2[i][m] = x1[i][m];
- else
- x2[i][m] = x1[i][m] + a*f[i][m];
- }
- }
+ x1 = s1->x;
+ x2 = s2->x;
- if (s2->flags & (1<<estCGP)) {
- /* Copy the CG p vector */
- x1 = s1->cg_p;
- x2 = s2->cg_p;
- for(i=start; i<end; i++)
- copy_rvec(x1[i],x2[i]);
- }
+#pragma omp parallel num_threads(gmx_omp_nthreads_get(emntUpdate))
+ {
+ int gf,i,m;
- if (DOMAINDECOMP(cr)) {
- s2->ddp_count = s1->ddp_count;
- if (s2->cg_gl_nalloc < s1->cg_gl_nalloc) {
- s2->cg_gl_nalloc = s1->cg_gl_nalloc;
- srenew(s2->cg_gl,s2->cg_gl_nalloc);
- }
- s2->ncg_gl = s1->ncg_gl;
- for(i=0; i<s2->ncg_gl; i++)
- s2->cg_gl[i] = s1->cg_gl[i];
- s2->ddp_count_cg_gl = s1->ddp_count_cg_gl;
- }
+ gf = 0;
+#pragma omp for schedule(static) nowait
+ for(i=start; i<end; i++)
+ {
+ if (md->cFREEZE)
+ {
+ gf = md->cFREEZE[i];
+ }
+ for(m=0; m<DIM; m++)
+ {
+ if (ir->opts.nFreeze[gf][m])
+ {
+ x2[i][m] = x1[i][m];
+ }
+ else
+ {
+ x2[i][m] = x1[i][m] + a*f[i][m];
+ }
+ }
+ }
- if (constr) {
- wallcycle_start(wcycle,ewcCONSTR);
- dvdlambda = 0;
- constrain(NULL,TRUE,TRUE,constr,&top->idef,
- ir,NULL,cr,count,0,md,
- s1->x,s2->x,NULL,s2->box,s2->lambda[efptBONDED],
- &dvdlambda,NULL,NULL,nrnb,econqCoord,FALSE,0,0);
- wallcycle_stop(wcycle,ewcCONSTR);
- }
+ if (s2->flags & (1<<estCGP))
+ {
+ /* Copy the CG p vector */
+ x1 = s1->cg_p;
+ x2 = s2->cg_p;
+#pragma omp for schedule(static) nowait
+ for(i=start; i<end; i++)
+ {
+ copy_rvec(x1[i],x2[i]);
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ s2->ddp_count = s1->ddp_count;
+ if (s2->cg_gl_nalloc < s1->cg_gl_nalloc)
+ {
+#pragma omp barrier
+ s2->cg_gl_nalloc = s1->cg_gl_nalloc;
+ srenew(s2->cg_gl,s2->cg_gl_nalloc);
+#pragma omp barrier
+ }
+ s2->ncg_gl = s1->ncg_gl;
+#pragma omp for schedule(static) nowait
+ for(i=0; i<s2->ncg_gl; i++)
+ {
+ s2->cg_gl[i] = s1->cg_gl[i];
+ }
+ s2->ddp_count_cg_gl = s1->ddp_count_cg_gl;
+ }
+ }
+
+ if (constr)
+ {
+ wallcycle_start(wcycle,ewcCONSTR);
+ dvdlambda = 0;
+ constrain(NULL,TRUE,TRUE,constr,&top->idef,
+ ir,NULL,cr,count,0,md,
+ s1->x,s2->x,NULL,bMolPBC,s2->box,
+ s2->lambda[efptBONDED],&dvdlambda,
+ NULL,NULL,nrnb,econqCoord,FALSE,0,0);
+ wallcycle_stop(wcycle,ewcCONSTR);
+ }
}
static void em_dd_partition_system(FILE *fplog,int step,t_commrec *cr,
ems->s.box,ems->s.x,&ems->s.hist,
ems->f,force_vir,mdatoms,enerd,fcd,
ems->s.lambda,graph,fr,vsite,mu_tot,t,NULL,NULL,TRUE,
- GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES | GMX_FORCE_VIRIAL |
+ GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
+ GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
(bNS ? GMX_FORCE_NS | GMX_FORCE_DOLR : 0));
/* Clear the unused shake virial and pressure */
dvdlambda = 0;
constrain(NULL,FALSE,FALSE,constr,&top->idef,
inputrec,NULL,cr,count,0,mdatoms,
- ems->s.x,ems->f,ems->f,ems->s.box,ems->s.lambda[efptBONDED],&dvdlambda,
+ ems->s.x,ems->f,ems->f,fr->bMolPBC,ems->s.box,
+ ems->s.lambda[efptBONDED],&dvdlambda,
NULL,&shake_vir,nrnb,econqForceDispl,FALSE,0,0);
if (fr->bSepDVDL && fplog)
fprintf(fplog,sepdvdlformat,"Constraints",t,dvdlambda);
}
/* Take a trial step (new coords in s_c) */
- do_em_step(cr,inputrec,mdatoms,s_min,c,s_min->s.cg_p,s_c,
- constr,top,nrnb,wcycle,-1);
+ do_em_step(cr,inputrec,mdatoms,fr->bMolPBC,s_min,c,s_min->s.cg_p,s_c,
+ constr,top,nrnb,wcycle,-1);
neval++;
/* Calculate energy for the trial step */
}
/* Take a trial step to this new point - new coords in s_b */
- do_em_step(cr,inputrec,mdatoms,s_min,b,s_min->s.cg_p,s_b,
- constr,top,nrnb,wcycle,-1);
+ do_em_step(cr,inputrec,mdatoms,fr->bMolPBC,s_min,b,s_min->s.cg_p,s_b,
+ constr,top,nrnb,wcycle,-1);
neval++;
/* Calculate energy for the trial step */
/* set new coordinates, except for first step */
if (count > 0) {
- do_em_step(cr,inputrec,mdatoms,s_min,stepsize,s_min->f,s_try,
- constr,top,nrnb,wcycle,count);
+ do_em_step(cr,inputrec,mdatoms,fr->bMolPBC,
+ s_min,stepsize,s_min->f,s_try,
+ constr,top,nrnb,wcycle,count);
}
evaluate_energy(fplog,bVerbose,cr,