From 9f544a7de5f29eeb62945cdaea6341b61177c827 Mon Sep 17 00:00:00 2001 From: Adam Herbst Date: Fri, 15 Apr 2011 10:27:30 -0400 Subject: [PATCH] A freeze group can now be allowed to move rigidly in some dimension by using "freezedim = R" in the .mdp file --- include/domdec.h | 6 ++++ include/types/inputrec.h | 10 ++++++ src/kernel/md.c | 5 ++- src/kernel/readir.c | 9 +++-- src/kernel/repl_ex.c | 30 +++++++++------- src/mdlib/domdec.c | 92 ++++++++++++++++++++++++++++++++++++++++++++++++ src/mdlib/pull.c | 14 ++++---- src/mdlib/pullutil.c | 4 ++- src/mdlib/sim_util.c | 62 ++++++++++++++++++++++++++++++++ src/mdlib/update.c | 17 ++++++--- 10 files changed, 221 insertions(+), 28 deletions(-) diff --git a/include/domdec.h b/include/domdec.h index 64b4d20a24..901ee4e8fc 100644 --- a/include/domdec.h +++ b/include/domdec.h @@ -141,6 +141,12 @@ void dd_atom_spread_real(gmx_domdec_t *dd,real v[]); void dd_atom_sum_real(gmx_domdec_t *dd,real v[]); /* Sum the contributions to a real for each atom over the neighboring cells. */ +void init_local_rigid_groups(t_inputrec *ir, t_mdatoms *mdatoms, t_state *state_local); +/* Initialize the rigid body groups */ + +void get_local_rigid_groups(t_inputrec *ir, t_mdatoms *mdatoms); +/* Update the rigid body groups according to the current domain decomposition */ + void dd_partition_system(FILE *fplog, gmx_large_int_t step, t_commrec *cr, diff --git a/include/types/inputrec.h b/include/types/inputrec.h index 7c993638fc..0cf1055124 100644 --- a/include/types/inputrec.h +++ b/include/types/inputrec.h @@ -146,6 +146,15 @@ typedef struct { } t_pull; typedef struct { + int *grpcnt; /* Number of home atoms in each freeze group */ + real *mass; /* Total mass of the home atoms of each freeze group */ + rvec *force; /* Total force on the home atoms of each freeze group for the last stepcnt steps */ + rvec *vel; /* Velocity of each freeze group, calculated using the force summed over all threads */ + double *dbuf; /* Buffer for summing over all threads */ + int stepcnt; /* Number of steps since the last summation and update */ +} t_rigid; + +typedef struct { int eI; /* Integration method */ gmx_large_int_t nsteps; /* number of steps to be taken */ int simulation_part; /* Used in checkpointing to separate chunks */ @@ -267,6 +276,7 @@ typedef struct { real wall_ewald_zfac; /* Scaling factor for the box for Ewald */ int ePull; /* Type of pulling: no, umbrella or constraint */ t_pull *pull; /* The data for center of mass pulling */ + t_rigid *rigid; /* The data for rigid body motion */ real cos_accel; /* Acceleration for viscosity calculation */ tensor deform; /* Triclinic deformation velocities (nm/ps) */ int userint1; /* User determined parameters */ diff --git a/src/kernel/md.c b/src/kernel/md.c index 1346ef9ad2..23a8e3367f 100644 --- a/src/kernel/md.c +++ b/src/kernel/md.c @@ -1367,7 +1367,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], dd_make_local_pull_groups(NULL,ir->pull,mdatoms); } } - + if (DOMAINDECOMP(cr)) { /* Distribute the charge groups over the nodes from the master node */ @@ -1379,6 +1379,9 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], } update_mdatoms(mdatoms,state->lambda); + + /* The rigid body groups can be initialized now that the atom data is there */ + init_local_rigid_groups(ir, mdatoms, state); if (MASTER(cr)) { diff --git a/src/kernel/readir.c b/src/kernel/readir.c index 1e241af77b..91d26be2d7 100644 --- a/src/kernel/readir.c +++ b/src/kernel/readir.c @@ -2026,10 +2026,15 @@ void do_index(const char* mdparin, const char *ndx, snew(ir->opts.nFreeze,nr); for(i=k=0; (iopts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0); + if(gmx_strncasecmp(ptr1[k],"Y",1)==0) //Y means the group is frozen in this dimension + ir->opts.nFreeze[i][j]=1; + else if(gmx_strncasecmp(ptr1[k],"R",1)==0) //R means the group is rigid in this dimension + ir->opts.nFreeze[i][j]=2; + else ir->opts.nFreeze[i][j]=0; + if (!ir->opts.nFreeze[i][j]) { if (gmx_strncasecmp(ptr1[k],"N",1) != 0) { - sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only " + sprintf(warnbuf,"Please use Y(ES) or N(O) or R(IGID) for freezedim only " "(not %s)", ptr1[k]); warning(wi,warn_buf); } diff --git a/src/kernel/repl_ex.c b/src/kernel/repl_ex.c index e4c9089951..8bff387895 100644 --- a/src/kernel/repl_ex.c +++ b/src/kernel/repl_ex.c @@ -89,20 +89,24 @@ static void repl_quantity(FILE *fplog,const gmx_multisim_t *ms, { bDiff = TRUE; } - } - if (bDiff) - { - if (re->type >= 0 && re->type < ereNR) - { - gmx_fatal(FARGS,"For replica exchange both %s and %s differ", - erename[re->type],erename[ere]); - } - /* Set the replica exchange type and quantities */ - re->type = ere; - snew(re->q,re->nrepl); - for(s=0; snsim; s++) + + if (bDiff) { - re->q[s] = qall[s]; + if (re->type >= 0 && re->type < ereNR) + { + gmx_fatal(FARGS,"For replica exchange both %s and %s differ", + erename[re->type],erename[ere]); + } + else + { + /* Set the replica exchange type and quantities */ + re->type = ere; + snew(re->q,re->nrepl); + for(s=0; snsim; s++) + { + re->q[s] = qall[s]; + } + } } } diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c index ca90e64d90..6fad81f005 100644 --- a/src/mdlib/domdec.c +++ b/src/mdlib/domdec.c @@ -8149,6 +8149,93 @@ void print_dd_statistics(t_commrec *cr,t_inputrec *ir,FILE *fplog) } } +void get_local_rigid_groups(t_inputrec *ir, t_mdatoms *mdatoms) { + /* redetermine the number and mass of local atoms in each rigid body group */ + int i, gf; + for(i = 0; i < ir->opts.ngfrz; i++) { + ir->rigid->grpcnt[i] = 0; + ir->rigid->mass[i] = 0; + } + for(i = 0; i < mdatoms->homenr; i++) { + gf = mdatoms->cFREEZE[i]; + ir->rigid->grpcnt[gf]++; + ir->rigid->mass[gf] += mdatoms->massT[i]; + } +} + +void init_local_rigid_groups(t_inputrec *ir, t_mdatoms *mdatoms, t_state *state_local) +{ + /* calculate the number of rigid body groups and total rigid dimensions */ + int ngrig = 0, ndrig = 0, rig, i, d; + for(i = 0; i < ir->opts.ngfrz; i++) { + rig = 0; + for(d = 0; d < DIM; d++) { + if(ir->opts.nFreeze[i][d] == 2) { + ndrig++; + rig = 1; + } + } + if(rig) ngrig++; + } + if(ngrig == 0) { + ir->rigid = NULL; + return; + } + + /* Initialize the rigid body state */ + snew(ir->rigid, 1); + t_rigid *rigid = ir->rigid; + snew(rigid->grpcnt, ir->opts.ngfrz); + snew(rigid->mass, ir->opts.ngfrz); + snew(rigid->force, ir->opts.ngfrz); + snew(rigid->vel, ir->opts.ngfrz); + rigid->stepcnt = 0; + + /* For summing, the buffer holds the mass and atom count of each rigid group, and force on each rigid dimension */ + snew(rigid->dbuf, ndrig + 2*ngrig); + + /* calculate the local mass of each group, initialize the force to zero, and set the velocity to the average of all atoms' initial velocities */ + for(i = 0; i < ir->opts.ngfrz; i++) { + rigid->grpcnt[i] = 0; + rigid->mass[i] = 0; + for(d = 0; d < DIM; d++) { + rigid->force[i][d] = 0; + rigid->vel[i][d] = 0; + } + } + for(i = 0; i < mdatoms->homenr; i++) { + int gf = mdatoms->cFREEZE[i]; + rigid->grpcnt[gf]++; + rigid->mass[gf] += mdatoms->massT[i]; + /* initialize the velocity of each rigid dimension to the average over all home atoms */ + // for(d = 0; d < DIM; d++) { + // if(ir->opts.nFreeze[gf][d] == 2) + // rigid->vel[gf][d] += state_local->v[i][d]; + // } + } + for(i = 0; i < ir->opts.ngfrz; i++) { + if(rigid->grpcnt[i] > 0) { + for(d = 0; d < DIM; d++) rigid->vel[i][d] /= rigid->grpcnt[i]; + } + } + /* char msg[500], app[500]; + sprintf(msg, "Thread %d has %d atoms and %d freeze groups:\n", cr->nodeid, mdatoms->homenr, ir->opts.ngfrz); + for(i = 0; i < ir->opts.ngfrz; i++) { + sprintf(app, "%d %d %d\tmass %f (%d)\tvel %f %f %f\tforce %f %f %f\n", + ir->opts.nFreeze[i][0], ir->opts.nFreeze[i][1], ir->opts.nFreeze[i][2], + ir->rigid->mass[i], ir->rigid->grpcnt[i], + ir->rigid->vel[i][0], ir->rigid->vel[i][1], ir->rigid->vel[i][2], + ir->rigid->force[i][0], ir->rigid->force[i][1], ir->rigid->force[i][2]); + strncat(msg, app, 500); + } + sprintf(app, "Address of grpcnt is %p\n", ir->rigid->grpcnt); + strncat(msg, app, 500); + printf("%s\n", msg); + */ + +} + + void dd_partition_system(FILE *fplog, gmx_large_int_t step, t_commrec *cr, @@ -8599,6 +8686,11 @@ void dd_partition_system(FILE *fplog, /* Update the local pull groups */ dd_make_local_pull_groups(dd,ir->pull,mdatoms); } + + if (ir->rigid) { + /* Update the local rigid body groups */ + get_local_rigid_groups(ir,mdatoms); + } add_dd_statistics(dd); diff --git a/src/mdlib/pull.c b/src/mdlib/pull.c index 548180094c..c6c14b7f95 100644 --- a/src/mdlib/pull.c +++ b/src/mdlib/pull.c @@ -1027,15 +1027,15 @@ static void make_local_pull_group(gmx_ga2la_t ga2la, if (ii >= start && ii < end) { /* This is a home atom, add it to the local pull group */ if (pg->nat_loc >= pg->nalloc_loc) { - pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1); - srenew(pg->ind_loc,pg->nalloc_loc); - if (pg->epgrppbc == epgrppbcCOS || pg->weight) { - srenew(pg->weight_loc,pg->nalloc_loc); - } + pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1); + srenew(pg->ind_loc,pg->nalloc_loc); + if (pg->epgrppbc == epgrppbcCOS || pg->weight) { + srenew(pg->weight_loc,pg->nalloc_loc); + } } pg->ind_loc[pg->nat_loc] = ii; if (pg->weight) { - pg->weight_loc[pg->nat_loc] = pg->weight[i]; + pg->weight_loc[pg->nat_loc] = pg->weight[i]; } pg->nat_loc++; } @@ -1183,7 +1183,7 @@ static void init_pull_group_index(FILE *fplog,t_commrec *cr, " For pulling the whole group will be frozen.\n\n", g); } - pg->invtm = 0.0; + pg->invtm = 1.0; pg->wscale = 1.0; } } diff --git a/src/mdlib/pullutil.c b/src/mdlib/pullutil.c index 38f42ef39f..d648b62b94 100644 --- a/src/mdlib/pullutil.c +++ b/src/mdlib/pullutil.c @@ -460,8 +460,10 @@ void pull_calc_coms(t_commrec *cr, g,wmass,wwmass,pgrp->invtm); } } +// printf("Pull group %d with %d atoms has wmass %f wwmass %f invtm %f, COM distance %f %f %f\n", +// g, pgrp->nat, wmass, wwmass, pgrp->invtm, pgrp->x[0], pgrp->x[1], pgrp->x[2]); } - + if (PULL_CYL(pull)) { /* Calculate the COMs for the cyclinder reference groups */ make_cyl_refgrps(cr,pull,md,pbc,t,x,xp); diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c index 60ef80e6c3..ffa5b81e43 100644 --- a/src/mdlib/sim_util.c +++ b/src/mdlib/sim_util.c @@ -822,6 +822,68 @@ void do_force(FILE *fplog,t_commrec *cr, } } + /* calculate the forces and increment the velocities of all rigid body groups */ +#ifndef RIGID_DEBUG + //#define RIGID_DEBUG +#endif + t_rigid *rigid = inputrec->rigid; + if (rigid && rigid->stepcnt > 0) { + /* add each rigid dimension to the summing buffer */ + int i,d,count = 0,hasRigid; + for(i = 0; i < inputrec->opts.ngfrz; i++) { + hasRigid = 0; + for(d = 0; d < DIM; d++) { + if(inputrec->opts.nFreeze[i][d] == 2) { + hasRigid = 1; + rigid->dbuf[count++] = rigid->force[i][d]; + } + } + if(hasRigid) { + rigid->dbuf[count++] = rigid->grpcnt[i]; + rigid->dbuf[count++] = rigid->mass[i]; +#ifdef RIGID_DEBUG + printf("%d (%d) BEFORE: group %d has force %f %f %f, vel %f %f %f, mass %f (%d)\n", + cr->nodeid, step, i, rigid->force[i][0], rigid->force[i][1], rigid->force[i][2], + rigid->vel[i][0], rigid->vel[i][1], rigid->vel[i][2], rigid->mass[i], rigid->grpcnt[i]); +#endif + } + } + /* sum them in parallel */ + if(cr && PAR(cr)) { + gmx_sumd(count, rigid->dbuf, cr); + } + /* extract the summed forces and masses, and increment the velocities */ + if(MASTER(cr)) { + int x = 3; + } + count = 0; + double grpMass, grpCnt; + for(i = 0; i < inputrec->opts.ngfrz; i++) { + hasRigid = 0; + for(d = 0; d < DIM; d++) { + if(inputrec->opts.nFreeze[i][d] == 2) { + hasRigid = 1; + rigid->force[i][d] = rigid->dbuf[count++]; + } + } + if(hasRigid) { + grpCnt = rigid->dbuf[count++]; + grpMass = rigid->dbuf[count++]; + for(d = 0; d < DIM; d++) { + /* damp the previous velocity to prevent unstable oscillations */ + rigid->vel[i][d] = 0.8 * rigid->vel[i][d] + rigid->force[i][d] / (grpMass * rigid->stepcnt); + } +#ifdef RIGID_DEBUG + printf("%d (%d) AFTER: group %d has force %f %f %f, vel %f %f %f, mass %f (%d)\n", + cr->nodeid, step, i, rigid->force[i][0], rigid->force[i][1], rigid->force[i][2], + rigid->vel[i][0], rigid->vel[i][1], rigid->vel[i][2], grpMass, (int)(grpCnt+.001)); +#endif + for(d = 0; d < DIM; d++) rigid->force[i][d] = 0; + } + } + rigid->stepcnt = 0; + } + if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F) { /* Calculate the center of mass forces, this requires communication, diff --git a/src/mdlib/update.c b/src/mdlib/update.c index 861b48766d..72ac7d4656 100644 --- a/src/mdlib/update.c +++ b/src/mdlib/update.c @@ -123,7 +123,7 @@ static void do_update_md(int start,int nrend,double dt, unsigned short ptype[],unsigned short cFREEZE[], unsigned short cACC[],unsigned short cTC[], rvec x[],rvec xprime[],rvec v[], - rvec f[],matrix M, + rvec f[],matrix M,t_rigid *rigid, gmx_bool bNH,gmx_bool bPR) { double imass,w_dt; @@ -216,11 +216,20 @@ static void do_update_md(int start,int nrend,double dt, } else { - v[n][d] = 0.0; - xprime[n][d] = x[n][d]; + /* if the atom is part of a rigid body in this dimension, + * use its group velocity and update its group force */ + if(nFreeze[gf][d] == 2) { + rigid->force[gf][d] += f[n][d]; + v[n][d] = rigid->vel[gf][d]; + xprime[n][d] = x[n][d] + v[n][d] * dt; + } else { /* otherwise do not move it at all */ + v[n][d] = 0.0; + xprime[n][d] = x[n][d]; + } } } } + if(rigid) rigid->stepcnt++; } } @@ -1648,7 +1657,7 @@ void update_coords(FILE *fplog, ekind->tcstat,ekind->grpstat,state->nosehoover_vxi, inputrec->opts.acc,inputrec->opts.nFreeze,md->invmass,md->ptype, md->cFREEZE,md->cACC,md->cTC, - state->x,xprime,state->v,force,M, + state->x,xprime,state->v,force,M,inputrec->rigid, bNH,bPR); } else -- 2.11.4.GIT