From 314440ab8282212da23c945d2c766d6b58a87b11 Mon Sep 17 00:00:00 2001 From: Carsten Kutzner Date: Tue, 23 Nov 2010 14:41:15 +0100 Subject: [PATCH] Enforced rotation: code cleanup --- include/pull_rotation.h | 4 +- include/types/inputrec.h | 4 +- src/gmxlib/tpxio.c | 2 +- src/gmxlib/txtdump.c | 2 +- src/kernel/readrot.c | 4 +- src/mdlib/domdec.c | 1 + src/mdlib/pull_rotation.c | 590 +++++++++++++++++++------------------- src/mdlib/sim_util.c | 2 +- 8 files changed, 300 insertions(+), 309 deletions(-) diff --git a/include/pull_rotation.h b/include/pull_rotation.h index 403d571c42..8de318ee22 100644 --- a/include/pull_rotation.h +++ b/include/pull_rotation.h @@ -109,7 +109,7 @@ extern void dd_make_local_rotation_groups(gmx_domdec_t *dd,t_rot *rot); * local arrays have to be updated (masses, shifts) */ extern void do_rotation(t_commrec *cr,t_inputrec *ir,matrix box,rvec x[],real t, - int step,gmx_wallcycle_t wcycle,gmx_bool bNS); + gmx_large_int_t step,gmx_wallcycle_t wcycle,gmx_bool bNS); /*! \brief Add the enforced rotation forces to the official force array. @@ -128,7 +128,7 @@ extern void do_rotation(t_commrec *cr,t_inputrec *ir,matrix box,rvec x[],real t, * \param step The time step, used for output. * \param t Time, used for output. */ -extern real add_rot_forces(t_rot *rot, rvec f[], t_commrec *cr, int step, real t); +extern real add_rot_forces(t_rot *rot, rvec f[], t_commrec *cr, gmx_large_int_t step, real t); /*! \brief Close the enforced rotation output files. diff --git a/include/types/inputrec.h b/include/types/inputrec.h index 7715b212bb..694fad55dd 100644 --- a/include/types/inputrec.h +++ b/include/types/inputrec.h @@ -171,8 +171,8 @@ typedef struct { typedef struct { int ngrp; /* Number of rotation groups */ - int nstrout; /* Main output frequency for angle and potential */ - int nsttout; /* Outfreq. for torque, fitangles, slab centers */ + int nstrout; /* Output frequency for main rotation outfile */ + int nstsout; /* Output frequency for per-slab data */ t_rotgrp *grp; /* Groups to rotate */ gmx_enfrot_t enfrot; /* Stores non-inputrec enforced rotation data */ } t_rot; diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c index a7b7054a8b..f220efb559 100644 --- a/src/gmxlib/tpxio.c +++ b/src/gmxlib/tpxio.c @@ -280,7 +280,7 @@ static void do_rot(t_fileio *fio, t_rot *rot,gmx_bool bRead, int file_version) gmx_fio_do_int(fio,rot->ngrp); gmx_fio_do_int(fio,rot->nstrout); - gmx_fio_do_int(fio,rot->nsttout); + gmx_fio_do_int(fio,rot->nstsout); if (bRead) snew(rot->grp,rot->ngrp); for(g=0; gngrp; g++) diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c index 4117af0654..8cc5c802c1 100644 --- a/src/gmxlib/txtdump.c +++ b/src/gmxlib/txtdump.c @@ -532,7 +532,7 @@ static void pr_rot(FILE *fp,int indent,t_rot *rot) int g; PI("rot_nstrout",rot->nstrout); - PI("rot_nsttout",rot->nsttout); + PI("rot_nstsout",rot->nstsout); PI("rot_ngrp",rot->ngrp); for(g=0; gngrp; g++) pr_rotgrp(fp,indent,g,&rot->grp[g]); diff --git a/src/kernel/readrot.c b/src/kernel/readrot.c index ff1c097cb4..66101513d5 100644 --- a/src/kernel/readrot.c +++ b/src/kernel/readrot.c @@ -71,8 +71,10 @@ extern char **read_rotparams(int *ninp_p,t_inpfile **inp_p,t_rot *rot, inp = *inp_p; /* read rotation parameters */ + CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group"); ITYPE("rot_nstrout", rot->nstrout, 100); - ITYPE("rot_nsttout", rot->nsttout, 1000); + CTYPE("Output frequency for per-slab data (angles, torques and slab centers)"); + ITYPE("rot_nstsout", rot->nstsout, 1000); CTYPE("Number of rotation groups"); ITYPE("rot_ngroups", rot->ngrp,1); diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c index 3544b55898..abc08743e2 100644 --- a/src/mdlib/domdec.c +++ b/src/mdlib/domdec.c @@ -41,6 +41,7 @@ #include "force.h" #include "pme.h" #include "pull.h" +#include "pull_rotation.h" #include "gmx_wallcycle.h" #include "mdrun.h" #include "nsgrid.h" diff --git a/src/mdlib/pull_rotation.c b/src/mdlib/pull_rotation.c index 16dab29526..6b19af91a9 100644 --- a/src/mdlib/pull_rotation.c +++ b/src/mdlib/pull_rotation.c @@ -111,7 +111,7 @@ typedef struct gmx_enfrot /* Global enforced rotation data for a single rotation group */ typedef struct gmx_enfrotgrp { - real degangle; /* Rotation angle in degree */ + real degangle; /* Rotation angle in degrees */ matrix rotmat; /* Rotation matrix */ atom_id *ind_loc; /* Local rotation indices */ int nat_loc; /* Number of local group atoms */ @@ -140,13 +140,15 @@ typedef struct gmx_enfrotgrp real *mc; /* Collective masses */ real *mc_sorted; real invmass; /* one over the total mass of the rotation group */ + + real torque_v; /* Torque in the direction of rotation vector */ + real angle_v; /* Actual angle of the whole rotation group */ /* Fixed rotation only */ + real weight_v; /* Weights for angle determination */ rvec *xr_loc; /* Local reference coords, correctly rotated */ rvec *x_loc_pbc; /* Local current coords, correct PBC image */ real *m_loc; /* Masses of the current local atoms */ - real fix_torque_v; /* Torque in the direction of rotation vector */ - real fix_angles_v; - real fix_weight_v; + /* Flexible rotation only */ int nslabs_alloc; /* For this many slabs memory is allocated */ int slab_first; /* Lowermost slab for that the calculation needs @@ -178,6 +180,22 @@ typedef struct gmx_enfrotgrp } t_gmx_enfrotgrp; +/* Activate output of forces for correctness checks */ +/* #define PRINT_FORCES */ +#ifdef PRINT_FORCES +#define PRINT_FORCE_J fprintf(stderr,"f%d = %15.8f %15.8f %15.8f\n",erg->xc_ref_ind[j],erg->f_rot_loc[j][XX], erg->f_rot_loc[j][YY], erg->f_rot_loc[j][ZZ]); +#define PRINT_POT_TAU fprintf(stderr,"potential = %15.8f\n" "torque = %15.8f", erg->V, erg->torque_v); +#else +#define PRINT_FORCE_J +#define PRINT_POT_TAU +#endif + +/* Shortcuts for often used queries */ +#define ISFLEX(rg) ( (rg->eType==erotgFLEX) || (rg->eType==erotgFLEXT) || (rg->eType==erotgFLEX2) || (rg->eType==erotgFLEX2T) ) +#define ISCOLL(rg) ( (rg->eType==erotgFLEX) || (rg->eType==erotgFLEXT) || (rg->eType==erotgFLEX2) || (rg->eType==erotgFLEX2T) || (rg->eType==erotgRMPF) || (rg->eType==erotgRM2PF) ) + + + static double** allocate_square_matrix(int dim) { int i; @@ -203,8 +221,8 @@ static void free_square_matrix(double** mat, int dim) } -/* Output rotation energy and torque for each rotation group */ -static void reduce_output(t_commrec *cr, t_rot *rot, real t) +/* Output rotation energy, torques, etc. for each rotation group */ +static void reduce_output(t_commrec *cr, t_rot *rot, real t, gmx_large_int_t step) { int g,i,islab,nslabs=0; int count; /* MPI element counter */ @@ -223,39 +241,23 @@ static void reduce_output(t_commrec *cr, t_rot *rot, real t) for (g=0; g < rot->ngrp; g++) { rotg = &rot->grp[g]; - erg=rotg->enfrotgrp; + erg = rotg->enfrotgrp; nslabs = erg->slab_last - erg->slab_first + 1; er->mpi_inbuf[count++] = erg->V; - switch (rotg->eType) + er->mpi_inbuf[count++] = erg->torque_v; + er->mpi_inbuf[count++] = erg->angle_v; + er->mpi_inbuf[count++] = erg->weight_v; /* weights are not needed for flex types, but this is just a single value */ + if (ISFLEX(rotg)) { - case erotgISO: - case erotgISOPF: - case erotgPM: - case erotgPMPF: - case erotgRM: - case erotgRMPF: - case erotgRM2: - case erotgRM2PF: - er->mpi_inbuf[count++] = erg->fix_torque_v; - er->mpi_inbuf[count++] = erg->fix_angles_v; - er->mpi_inbuf[count++] = erg->fix_weight_v; - break; - case erotgFLEX: - case erotgFLEXT: - case erotgFLEX2: - case erotgFLEX2T: - /* (Re-)allocate memory for MPI buffer: */ - if (er->mpi_bufsize < count+nslabs) - { - er->mpi_bufsize = count+nslabs; - srenew(er->mpi_inbuf , er->mpi_bufsize); - srenew(er->mpi_outbuf, er->mpi_bufsize); - } - for (i=0; impi_inbuf[count++] = erg->slab_torque_v[i]; - break; - default: - break; + /* (Re-)allocate memory for MPI buffer: */ + if (er->mpi_bufsize < count+nslabs) + { + er->mpi_bufsize = count+nslabs; + srenew(er->mpi_inbuf , er->mpi_bufsize); + srenew(er->mpi_outbuf, er->mpi_bufsize); + } + for (i=0; impi_inbuf[count++] = erg->slab_torque_v[i]; } } #ifdef GMX_MPI @@ -268,32 +270,16 @@ static void reduce_output(t_commrec *cr, t_rot *rot, real t) for (g=0; g < rot->ngrp; g++) { rotg = &rot->grp[g]; - erg=rotg->enfrotgrp; + erg = rotg->enfrotgrp; nslabs = erg->slab_last - erg->slab_first + 1; - erg->V = er->mpi_outbuf[count++]; - switch (rotg->eType) + erg->V = er->mpi_outbuf[count++]; + erg->torque_v = er->mpi_outbuf[count++]; + erg->angle_v = er->mpi_outbuf[count++]; + erg->weight_v = er->mpi_outbuf[count++]; + if (ISFLEX(rotg)) { - case erotgISO: - case erotgISOPF: - case erotgPM: - case erotgPMPF: - case erotgRM: - case erotgRMPF: - case erotgRM2: - case erotgRM2PF: - erg->fix_torque_v = er->mpi_outbuf[count++]; - erg->fix_angles_v = er->mpi_outbuf[count++]; - erg->fix_weight_v = er->mpi_outbuf[count++]; - break; - case erotgFLEX: - case erotgFLEXT: - case erotgFLEX2: - case erotgFLEX2T: - for (i=0; islab_torque_v[i] = er->mpi_outbuf[count++]; - break; - default: - break; + for (i=0; islab_torque_v[i] = er->mpi_outbuf[count++]; } } } @@ -302,26 +288,27 @@ static void reduce_output(t_commrec *cr, t_rot *rot, real t) /* Output */ if (MASTER(cr)) { - /* Av. angle and total torque for each rotation group */ + /* Angle and torque for each rotation group */ for (g=0; g < rot->ngrp; g++) { rotg=&rot->grp[g]; - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); + bFlex = ISFLEX(rotg); erg=rotg->enfrotgrp; - /* Output to main rotation log file: */ - if (!bFlex) + /* Output to main rotation output file: */ + if ( do_per_step(step, rot->nstrout) ) { - fprintf(er->out_rot, "%12.4f%12.3e", - (erg->fix_angles_v/erg->fix_weight_v)*180.0*M_1_PI, - erg->fix_torque_v); + if (bFlex) + fprintf(er->out_rot, "%12.4f", erg->angle_v); /* RMSD fit angle */ + else + fprintf(er->out_rot, "%12.4f", (erg->angle_v/erg->weight_v)*180.0*M_1_PI); + fprintf(er->out_rot, "%12.3e", erg->torque_v); + fprintf(er->out_rot, "%12.3e", erg->V); } - fprintf(er->out_rot, "%12.3e", erg->V); - + /* Output to torque log file: */ - if (bFlex) + if ( bFlex && do_per_step(step, rot->nstsout) ) { fprintf(er->out_torque, "%12.3e%6d", t, g); for (i=erg->slab_first; i<=erg->slab_last; i++) @@ -334,14 +321,15 @@ static void reduce_output(t_commrec *cr, t_rot *rot, real t) fprintf(er->out_torque , "\n"); } } - fprintf(er->out_rot, "\n"); + if ( do_per_step(step, rot->nstrout) ) + fprintf(er->out_rot, "\n"); } } /* Add the forces from enforced rotation potential to the local forces. * Should be called after the SR forces have been evaluated */ -extern real add_rot_forces(t_rot *rot, rvec f[], t_commrec *cr, int step, real t) +extern real add_rot_forces(t_rot *rot, rvec f[], t_commrec *cr, gmx_large_int_t step, real t) { int g,l,ii; t_rotgrp *rotg; @@ -355,8 +343,8 @@ extern real add_rot_forces(t_rot *rot, rvec f[], t_commrec *cr, int step, real t /* Reduce energy,torque, angles etc. to get the sum values (per rotation group) * on the master and output these values to file. */ - if (do_per_step(step, rot->nsttout)) - reduce_output(cr, rot, t); + if ( do_per_step(step, rot->nstrout) || do_per_step(step, rot->nstsout) ) + reduce_output(cr, rot, t, step); /* Total rotation potential is the sum over all rotation groups */ er->Vrot = 0.0; @@ -612,17 +600,6 @@ static void print_aligned_short(FILE *fp, char *str) } -/* Right-aligned output of value with standard width */ -static void print_aligned_group(FILE *fp, char *str, int g) -{ - char sbuf[STRLEN]; - - - sprintf(sbuf, "%s%d", str, g); - fprintf(fp, "%12s", sbuf); -} - - static FILE *open_output_file(const char *fn, int steps, const char what[]) { FILE *fp; @@ -630,7 +607,7 @@ static FILE *open_output_file(const char *fn, int steps, const char what[]) fp = ffopen(fn, "w"); - fprintf(fp, "# Output of %s is written at intervals of %d time step%s.\n#\n", + fprintf(fp, "# Output of %s is written in intervals of %d time step%s.\n#\n", what,steps, steps>1 ? "s":""); return fp; @@ -648,11 +625,10 @@ static FILE *open_slab_out(t_rot *rot, const char *fn) for (g=0; gngrp; g++) { rotg = &rot->grp[g]; - if ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ) + if (ISFLEX(rotg)) { if (NULL == fp) - fp = open_output_file(fn, rot->nsttout, "gaussian weighted slab centers"); + fp = open_output_file(fn, rot->nstsout, "gaussian weighted slab centers"); fprintf(fp, "# Rotation group %d (%s), slab distance %f nm, %s.\n", g, erotg_names[rotg->eType], rotg->slab_dist, rotg->bMassW? "centers of mass":"geometrical centers"); } @@ -690,7 +666,7 @@ static FILE *open_rot_out(const char *fn, t_rot *rot, const output_env_t oenv, int g,nsets; t_rotgrp *rotg; const char **setname; - char buf[50]; + char buf[50], buf2[75]; gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ gmx_bool bFlex; @@ -701,22 +677,24 @@ static FILE *open_rot_out(const char *fn, t_rot *rot, const output_env_t oenv, } else { - fp = xvgropen(fn, "Rotation angles and energy", "Time [ps]", "angles [degree] and energies [kJ/mol]", oenv); + fp = xvgropen(fn, "Rotation angles and energy", "Time [ps]", "angles [degrees] and energies [kJ/mol]", oenv); + fprintf(fp, "# Output of enforced rotation data is written in intervals of %d time step%s.\n#\n", rot->nstrout, rot->nstrout > 1 ? "s":""); fprintf(fp, "# The scalar tau is the torque [kJ/mol] in the direction of the rotation vector v.\n"); - fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot_vec.\n#\n"); + fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot_vec.\n"); + fprintf(fp, "# For flexible groups, tau(t,n) from all slabs n have been summed in a single value tau(t) here.\n"); + fprintf(fp, "# The torques tau(t,n) are found in the rottorque.log (-rt) output file\n"); + fprintf(fp, "#\n"); for (g=0; gngrp; g++) { rotg = &rot->grp[g]; erg=rotg->enfrotgrp; - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); - + bFlex = ISFLEX(rotg); fprintf(fp, "# Rotation group %d, potential type '%s':\n" , g, erotg_names[rotg->eType]); fprintf(fp, "# rot_massw%d %s\n" , g, yesno_names[rotg->bMassW]); fprintf(fp, "# rot_vec%d %12.5e %12.5e %12.5e\n" , g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]); - fprintf(fp, "# rot_rate%d %12.5e degree/ps\n" , g, rotg->rate); + fprintf(fp, "# rot_rate%d %12.5e degrees/ps\n" , g, rotg->rate); fprintf(fp, "# rot_k%d %12.5e kJ/(mol*nm^2)\n" , g, rotg->k); if ( rotg->eType==erotgISO || rotg->eType==erotgPM || rotg->eType==erotgRM || rotg->eType==erotgRM2) fprintf(fp, "# rot_pivot%d %12.5e %12.5e %12.5e nm\n", g, rotg->pivot[XX], rotg->pivot[YY], rotg->pivot[ZZ]); @@ -755,33 +733,38 @@ static FILE *open_rot_out(const char *fn, t_rot *rot, const output_env_t oenv, for (g=0; gngrp; g++) { rotg = &rot->grp[g]; - sprintf(buf, "theta_ref%d [degree]", g); - print_aligned_group(fp, "theta_ref", g); - setname[nsets] = strdup(buf); + sprintf(buf, "theta_ref%d", g); + print_aligned(fp, buf); + sprintf(buf2, "%s [degrees]", buf); + setname[nsets] = strdup(buf2); nsets++; } for (g=0; gngrp; g++) { rotg = &rot->grp[g]; - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); + bFlex = ISFLEX(rotg); /* For flexible axis rotation we use RMSD fitting to determine the * actual angle of the rotation group */ - if (!bFlex) - { - sprintf(buf, "theta-av%d [degree]", g); - print_aligned_group(fp, "theta_av", g); - setname[nsets] = strdup(buf); - nsets++; - sprintf(buf, "tau%d [kJ/mol]", g); - print_aligned_group(fp, "tau", g); - setname[nsets] = strdup(buf); - nsets++; - } - sprintf(buf, "energy%d [kJ/mol]", g); - print_aligned_group(fp, "energy", g); - setname[nsets] = strdup(buf); + if (bFlex) + sprintf(buf, "theta-fit%d", g); + else + sprintf(buf, "theta-av%d", g); + print_aligned(fp, buf); + sprintf(buf2, "%s [degrees]", buf); + setname[nsets] = strdup(buf2); + nsets++; + + sprintf(buf, "tau%d", g); + print_aligned(fp, buf); + sprintf(buf2, "%s [kJ/mol]", buf); + setname[nsets] = strdup(buf2); + nsets++; + + sprintf(buf, "energy%d", g); + print_aligned(fp, buf); + sprintf(buf2, "%s [kJ/mol]", buf); + setname[nsets] = strdup(buf2); nsets++; } fprintf(fp, "\n#\n"); @@ -806,13 +789,12 @@ static FILE *open_angles_out(t_rot *rot, const char *fn) /* Open output file and write some information about it's structure: */ - fp = open_output_file(fn, rot->nstrout, "rotation group angles"); + fp = open_output_file(fn, rot->nstsout, "rotation group angles"); fprintf(fp, "# All angles given in degrees, time in ps.\n"); for (g=0; gngrp; g++) { rotg = &rot->grp[g]; - if ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ) + if (ISFLEX(rotg)) { fprintf(fp, "# Rotation group %d (%s), slab distance %f nm, fit type %s.\n", g, erotg_names[rotg->eType], rotg->slab_dist, erotg_fitnames[rotg->eFittype]); @@ -823,7 +805,6 @@ static FILE *open_angles_out(t_rot *rot, const char *fn) print_aligned_short(fp, "t"); print_aligned_short(fp, "grp"); print_aligned(fp, "theta_ref"); - print_aligned(fp, "theta_fit"); print_aligned_short(fp, "slab"); print_aligned_short(fp, "atoms"); print_aligned(fp, "theta_fit"); @@ -845,13 +826,12 @@ static FILE *open_torque_out(t_rot *rot, const char *fn) t_rotgrp *rotg; - fp = open_output_file(fn, rot->nsttout,"torques"); + fp = open_output_file(fn, rot->nstsout,"torques"); for (g=0; gngrp; g++) { rotg = &rot->grp[g]; - if ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ) + if (ISFLEX(rotg)) { fprintf(fp, "# Rotation group %d (%s), slab distance %f nm.\n", g, erotg_names[rotg->eType], rotg->slab_dist); fprintf(fp, "# The scalar tau is the torque [kJ/mol] in the direction of the rotation vector.\n"); @@ -1013,7 +993,7 @@ static void weigh_coords(rvec* str, real* weight, int natoms) } -static double opt_angle_analytic( +static real opt_angle_analytic( rvec* ref_s, rvec* act_s, real* weight, @@ -1166,13 +1146,61 @@ static double opt_angle_analytic( sfree(eigvec[i]); sfree(eigvec); - return opt_angle; + return (real) opt_angle; } -/* Determine actual angle of this slab by RMSD fit to the reference */ +/* Determine angle of the group by RMSD fit to the reference */ /* Not parallelized, call this routine only on the master */ -static void flex_fit_angle( +static real flex_fit_angle(t_rotgrp *rotg) +{ + int i; + rvec *fitcoords=NULL; + rvec center; /* Center of positions passed to the fit routine */ + real fitangle; /* Angle of the rotation group derived by fitting */ + rvec coord; + real scal; + gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ + + + erg=rotg->enfrotgrp; + + /* Get the center of the rotation group. + * Note, again, erg->xc has been sorted in do_flexible */ + get_center(erg->xc, erg->mc_sorted, rotg->nat, center); + + /* === Determine the optimal fit angle for the rotation group === */ + if (rotg->eFittype == erotgFitNORM) + { + /* Normalize every position to it's reference length */ + for (i=0; inat; i++) + { + /* Put the center of the positions into the origin */ + rvec_sub(erg->xc[i], center, coord); + /* Determine the scaling factor for the length: */ + scal = erg->xc_ref_length[erg->xc_sortind[i]] / norm(coord); + /* Get position, multiply with the scaling factor and save */ + svmul(scal, coord, erg->xc_norm[i]); + } + fitcoords = erg->xc_norm; + } + else + { + fitcoords = erg->xc; + } + /* From the point of view of the current positions, the reference has rotated + * backwards. Since we output the angle relative to the fixed reference, + * we need the minus sign. */ + fitangle = -opt_angle_analytic(erg->xc_ref_sorted, fitcoords, erg->mc_sorted, + rotg->nat, erg->xc_ref_center, center, rotg->vec); + + return fitangle; +} + + +/* Determine actual angle of each slab by RMSD fit to the reference */ +/* Not parallelized, call this routine only on the master */ +static void flex_fit_angle_perslab( int g, t_rotgrp *rotg, double t, @@ -1181,19 +1209,16 @@ static void flex_fit_angle( { int i,l,n,islab,ind; rvec curr_x, ref_x; - rvec *fitcoords=NULL; rvec act_center; /* Center of actual positions that are passed to the fit routine */ rvec ref_center; /* Same for the reference positions */ - double fitangle; /* This will be the actual angle of the rotation group derived - * from an RMSD fit to the reference structure at t=0 */ + real fitangle; /* Angle of a slab derived from an RMSD fit to + * the reference structure at t=0 */ t_gmx_slabdata *sd; - rvec coord; - real scal; gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ real OOm_av; /* 1/average_mass of a rotation group atom */ real m_rel; /* Relative mass of a rotation group atom */ - + erg=rotg->enfrotgrp; /* Average mass of a rotation group atom: */ @@ -1210,23 +1235,23 @@ static void flex_fit_angle( sd = &(rotg->enfrotgrp->slab_data[islab]); sd->nat = erg->lastatom[islab]-erg->firstatom[islab]+1; ind = 0; - + /* Loop over the relevant atoms in the slab */ for (l=erg->firstatom[islab]; l<=erg->lastatom[islab]; l++) - { + { /* Current position of this atom: x[ii][XX/YY/ZZ] */ copy_rvec(erg->xc[l], curr_x); - + /* The (unrotated) reference position of this atom is copied to ref_x. * Beware, the xc coords have been sorted in do_flexible */ copy_rvec(erg->xc_ref_sorted[l], ref_x); - + /* Save data for doing angular RMSD fit later */ /* Save the current atom position */ copy_rvec(curr_x, sd->x[ind]); /* Save the corresponding reference position */ copy_rvec(ref_x , sd->ref[ind]); - + /* Maybe also mass-weighting was requested. If yes, additionally * multiply the weights with the relative mass of the atom. If not, * multiply with unity. */ @@ -1240,40 +1265,11 @@ static void flex_fit_angle( } } - /* Get the center of the whole rotation group. Note, again, the erg->xc have - * been sorted in do_flexible */ - get_center(erg->xc, erg->mc_sorted, rotg->nat, act_center); - /******************************/ /* Now do the fit calculation */ /******************************/ - /* === Determine the optimal fit angle for the whole rotation group === */ - if (rotg->eFittype == erotgFitNORM) - { - /* Normalize every position to it's reference length - * prior to performing the fit */ - for (i=0; inat; i++) - { - /* First put the center of the positions into the origin */ - rvec_sub(erg->xc[i], act_center, coord); - /* Determine the scaling factor for the length: */ - scal = erg->xc_ref_length[erg->xc_sortind[i]] / norm(coord); - /* Get position, multiply with the scaling factor and save in buf[i] */ - svmul(scal, coord, erg->xc_norm[i]); - } - fitcoords = erg->xc_norm; - } - else - { - fitcoords = erg->xc; - } - /* Note that from the point of view of the current positions, the reference has rotated backwards, - * but we want to output the angle relative to the fixed reference, therefore the minus sign. */ - fitangle = -opt_angle_analytic(erg->xc_ref_sorted, fitcoords, erg->mc_sorted, - rotg->nat, erg->xc_ref_center, act_center, rotg->vec); - fprintf(fp, "%12.3e%6d%12.3f%12.3lf", t, g, degangle, fitangle); - + fprintf(fp, "%12.3e%6d%12.3f", t, g, degangle); /* === Now do RMSD fitting for each slab === */ /* We require at least SLAB_MIN_ATOMS in a slab, such that the fit makes sense. */ @@ -1290,7 +1286,7 @@ static void flex_fit_angle( get_center(sd->x , sd->weight, sd->nat, act_center); if (rotg->eFittype == erotgFitNORM) { - /* Normalize every position to it's reference length + /* Normalize every position to it's reference length * prior to performing the fit */ for (i=0; inat;i++) /* Center */ { @@ -1303,7 +1299,8 @@ static void flex_fit_angle( clear_rvec(ref_center); clear_rvec(act_center); } - fitangle = -opt_angle_analytic(sd->ref, sd->x, sd->weight, sd->nat, ref_center, act_center, rotg->vec); + fitangle = -opt_angle_analytic(sd->ref, sd->x, sd->weight, sd->nat, + ref_center, act_center, rotg->vec); fprintf(fp, "%6d%6d%12.3f", n, sd->nat, fitangle); } } @@ -1767,22 +1764,16 @@ static real do_flex2_lowlevel( for (m=0; mf_rot_loc[j][m] = rotg->k * (-sum1vec[m] + sum2vec[m] - sum3vec[m] + 0.5*sum4vec[m]); -#ifdef INFOF - fprintf(stderr," FORCE on ATOM %d/%d = (%15.8f %15.8f %15.8f) \n", - j,ii,erg->f_rot_loc[j][XX], erg->f_rot_loc[j][YY], erg->f_rot_loc[j][ZZ]); -#endif - #ifdef SUM_PARTS fprintf(stderr, "sum1: %15.8f %15.8f %15.8f\n", -rotg->k*sum1vec[XX], -rotg->k*sum1vec[YY], -rotg->k*sum1vec[ZZ]); fprintf(stderr, "sum2: %15.8f %15.8f %15.8f\n", rotg->k*sum2vec[XX], rotg->k*sum2vec[YY], rotg->k*sum2vec[ZZ]); fprintf(stderr, "sum3: %15.8f %15.8f %15.8f\n", -rotg->k*sum3vec[XX], -rotg->k*sum3vec[YY], -rotg->k*sum3vec[ZZ]); fprintf(stderr, "sum4: %15.8f %15.8f %15.8f\n", 0.5*rotg->k*sum4vec[XX], 0.5*rotg->k*sum4vec[YY], 0.5*rotg->k*sum4vec[ZZ]); #endif - } /* END of loop over local atoms */ -#ifdef INFOF - fprintf(stderr, "THE POTENTIAL IS V=%f\n", V); -#endif + PRINT_FORCE_J + + } /* END of loop over local atoms */ return V; } @@ -1948,12 +1939,9 @@ static real do_flex_lowlevel( * array after the normal forces have been evaluated */ for(m=0; mf_rot_loc[j][m] = rotg->k*tmp_f[m]; -#ifdef INFOF - fprintf(stderr," FORCE on atom %d = %15.8f %15.8f %15.8f 1: %15.8f %15.8f %15.8f 2: %15.8f %15.8f %15.8f\n", iigrp, - rotg->k*tmp_f[XX] , rotg->k*tmp_f[YY] , rotg->k*tmp_f[ZZ] , - -rotg->k*sum_n1[XX], -rotg->k*sum_n1[YY], -rotg->k*sum_n1[ZZ], - rotg->k*sum_n2[XX], rotg->k*sum_n2[YY], rotg->k*sum_n2[ZZ]); -#endif + + PRINT_FORCE_J + } /* END of loop over local atoms */ return V; @@ -2152,14 +2140,15 @@ static void get_firstlast_slab_check( /* Enforced rotation with a flexible axis */ static void do_flexible( t_commrec *cr, - gmx_enfrot_t enfrot, /* Other rotation data */ - t_rotgrp *rotg, /* The rotation group */ - int g, /* Group number */ - rvec x[], /* The local positions */ + gmx_enfrot_t enfrot, /* Other rotation data */ + t_rotgrp *rotg, /* The rotation group */ + int g, /* Group number */ + rvec x[], /* The local positions */ matrix box, - double t, /* Time in picoseconds */ - int step, /* The time step */ - gmx_bool bOutstep) + double t, /* Time in picoseconds */ + gmx_large_int_t step, /* The time step */ + gmx_bool bOutstepRot, /* Output to main rotation output file */ + gmx_bool bOutstepSlab) /* Output per-slab data */ { int l,nslabs; real sigma; /* The Gaussian width sigma */ @@ -2172,8 +2161,7 @@ static void do_flexible( sigma = 0.7*rotg->slab_dist; /* Sort the collective coordinates erg->xc along the rotation vector. This is - * an optimization for the inner loop. - */ + * an optimization for the inner loop. */ sort_collective_coordinates(rotg, enfrot->data); /* Determine the first relevant slab for the first atom and the last @@ -2185,7 +2173,7 @@ static void do_flexible( get_firstlast_atom_per_slab(rotg, cr); /* Determine the gaussian-weighted center of positions for all slabs */ - get_slab_centers(rotg,erg->xc,erg->mc_sorted,cr,g,t,enfrot->out_slabs,bOutstep,FALSE); + get_slab_centers(rotg,erg->xc,erg->mc_sorted,cr,g,t,enfrot->out_slabs,bOutstepSlab,FALSE); /* Clear the torque per slab from last time step: */ nslabs = erg->slab_last - erg->slab_first + 1; @@ -2195,17 +2183,35 @@ static void do_flexible( /* Call the rotational forces kernel */ GMX_MPE_LOG(ev_flexll_start); if (rotg->eType == erotgFLEX || rotg->eType == erotgFLEXT) - erg->V = do_flex_lowlevel(rotg, sigma, x, bOutstep, box, cr); + erg->V = do_flex_lowlevel(rotg, sigma, x, bOutstepRot, box, cr); else if (rotg->eType == erotgFLEX2 || rotg->eType == erotgFLEX2T) - erg->V = do_flex2_lowlevel(rotg, sigma, x, bOutstep, box, cr); + erg->V = do_flex2_lowlevel(rotg, sigma, x, bOutstepRot, box, cr); else gmx_fatal(FARGS, "Unknown flexible rotation type"); GMX_MPE_LOG(ev_flexll_finish); - /* Determine actual angle of this slab by RMSD fit and output to file - Let's hope */ - /* this only happens once in a while, since this is not parallelized! */ - if (bOutstep && MASTER(cr)) - flex_fit_angle(g, rotg, t, erg->degangle, enfrot->out_angles); + /* Determine angle by RMSD fit to the reference - Let's hope this */ + /* only happens once in a while, since this is not parallelized! */ + if (MASTER(cr)) + { + if (bOutstepRot) + { + /* Fit angle of the whole rotation group */ + erg->angle_v = flex_fit_angle(rotg); + } + if (bOutstepSlab) + { + /* Fit angle of each slab */ + flex_fit_angle_perslab(g, rotg, t, erg->degangle, enfrot->out_angles); + } + } + + /* Lump together the torques from all slabs: */ + erg->torque_v = 0.0; + for (l=0; ltorque_v += erg->slab_torque_v[l]; + + PRINT_POT_TAU } @@ -2266,10 +2272,10 @@ static void do_fixed( rvec x[], /* The positions */ matrix box, /* The simulation box */ double t, /* Time in picoseconds */ - int step, /* The time step */ + gmx_large_int_t step, /* The time step */ gmx_bool bTorque) { - int i,m; + int j,m; rvec dr; rvec tmp_f; /* Force */ real alpha; /* a single angle between an actual and a reference position */ @@ -2288,28 +2294,22 @@ static void do_fixed( erg=rotg->enfrotgrp; bProject = (rotg->eType==erotgPM) || (rotg->eType==erotgPMPF); - /* Clear values from last time step */ - erg->V = 0.0; - erg->fix_torque_v = 0.0; - erg->fix_angles_v = 0.0; - erg->fix_weight_v = 0.0; - N_M = rotg->nat * erg->invmass; /* Each process calculates the forces on its local atoms */ - for (i=0; inat_loc; i++) + for (j=0; jnat_loc; j++) { /* Calculate (x_i-x_c) resp. (x_i-u) */ - rvec_sub(erg->x_loc_pbc[i], erg->xc_center, tmpvec); + rvec_sub(erg->x_loc_pbc[j], erg->xc_center, tmpvec); /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */ - rvec_sub(erg->xr_loc[i], tmpvec, dr); + rvec_sub(erg->xr_loc[j], tmpvec, dr); if (bProject) project_onto_plane(dr, rotg->vec); /* Mass-weighting */ - wi = N_M*erg->m_loc[i]; + wi = N_M*erg->m_loc[j]; /* Store the additional force so that it can be added to the force * array after the normal forces have been evaluated */ @@ -2317,19 +2317,19 @@ static void do_fixed( for (m=0; mf_rot_loc[i][m] = tmp_f[m]; + erg->f_rot_loc[j][m] = tmp_f[m]; erg->V += 0.5*k_wi*sqr(dr[m]); } if (bTorque) { /* Add to the torque of this rotation group */ - erg->fix_torque_v += torque(rotg->vec, tmp_f, erg->x_loc_pbc[i], erg->xc_center); + erg->torque_v += torque(rotg->vec, tmp_f, erg->x_loc_pbc[j], erg->xc_center); /* Calculate the angle between reference and actual rotation group atom. */ - angle(rotg, tmpvec, erg->xr_loc[i], &alpha, &weight); /* angle in rad, weighted */ - erg->fix_angles_v += alpha * weight; - erg->fix_weight_v += weight; + angle(rotg, tmpvec, erg->xr_loc[j], &alpha, &weight); /* angle in rad, weighted */ + erg->angle_v += alpha * weight; + erg->weight_v += weight; } /* If you want enforced rotation to contribute to the virial, * activate the following lines: @@ -2341,11 +2341,12 @@ static void do_fixed( vir[j][m] += 0.5*f[ii][j]*dr[m]; } */ -#ifdef INFOF - fprintf(stderr,"step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step, cr->nodeid, - erg->xc_ref_ind[i],erg->f_rot_loc[i][XX], erg->f_rot_loc[i][YY], erg->f_rot_loc[i][ZZ],erg->fix_torque_v); -#endif + + PRINT_FORCE_J + } /* end of loop over local rotation group atoms */ + + PRINT_POT_TAU } @@ -2356,7 +2357,7 @@ static void do_radial_motion( rvec x[], /* The positions */ matrix box, /* The simulation box */ double t, /* Time in picoseconds */ - int step, /* The time step */ + gmx_large_int_t step, /* The time step */ gmx_bool bTorque) { int j; @@ -2366,7 +2367,7 @@ static void do_radial_motion( gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ rvec xj_u; /* xj - u */ rvec tmpvec; - real fac,fac2,sum; + real fac,fac2,sum=0.0; rvec pj; /* For mass weighting: */ @@ -2376,13 +2377,6 @@ static void do_radial_motion( erg=rotg->enfrotgrp; - /* Clear values from last time step */ - erg->V = 0.0; - sum = 0.0; - erg->fix_torque_v = 0.0; - erg->fix_angles_v = 0.0; - erg->fix_weight_v = 0.0; - N_M = rotg->nat * erg->invmass; /* Each process calculates the forces on its local atoms */ @@ -2412,19 +2406,21 @@ static void do_radial_motion( if (bTorque) { /* Add to the torque of this rotation group */ - erg->fix_torque_v += torque(rotg->vec, tmp_f, erg->x_loc_pbc[j], erg->xc_center); + erg->torque_v += torque(rotg->vec, tmp_f, erg->x_loc_pbc[j], erg->xc_center); /* Calculate the angle between reference and actual rotation group atom. */ angle(rotg, xj_u, erg->xr_loc[j], &alpha, &weight); /* angle in rad, weighted */ - erg->fix_angles_v += alpha * weight; - erg->fix_weight_v += weight; + erg->angle_v += alpha * weight; + erg->weight_v += weight; } -#ifdef INFOF - fprintf(stderr,"RM: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step, cr->nodeid, - erg->xc_ref_ind[j],erg->f_rot_loc[j][XX], erg->f_rot_loc[j][YY], erg->f_rot_loc[j][ZZ],erg->fix_torque_v); -#endif + + PRINT_FORCE_J + } /* end of loop over local rotation group atoms */ erg->V = 0.5*rotg->k*sum; + + PRINT_POT_TAU + } @@ -2435,7 +2431,7 @@ static void do_radial_motion_pf( rvec x[], /* The positions */ matrix box, /* The simulation box */ double t, /* Time in picoseconds */ - int step, /* The time step */ + gmx_large_int_t step, /* The time step */ gmx_bool bTorque) { int i,ii,iigrp,j; @@ -2449,7 +2445,7 @@ static void do_radial_motion_pf( rvec tmpvec, tmpvec2; rvec innersumvec; /* Precalculation of the inner sum */ rvec innersumveckM; - real fac,fac2,V; + real fac,fac2,V=0.0; rvec qi,qj; /* For mass weighting: */ @@ -2459,13 +2455,6 @@ static void do_radial_motion_pf( erg=rotg->enfrotgrp; - /* Clear values from last time step */ - erg->V = 0.0; - V = 0.0; - erg->fix_torque_v = 0.0; - erg->fix_angles_v = 0.0; - erg->fix_weight_v = 0.0; - N_M = rotg->nat * erg->invmass; /* Get the current center of the rotation group: */ @@ -2542,19 +2531,20 @@ static void do_radial_motion_pf( if (bTorque) { /* Add to the torque of this rotation group */ - erg->fix_torque_v += torque(rotg->vec, tmp_f, xj, erg->xc_center); + erg->torque_v += torque(rotg->vec, tmp_f, xj, erg->xc_center); /* Calculate the angle between reference and actual rotation group atom. */ angle(rotg, xj_xc, yj0_yc0, &alpha, &weight); /* angle in rad, weighted */ - erg->fix_angles_v += alpha * weight; - erg->fix_weight_v += weight; + erg->angle_v += alpha * weight; + erg->weight_v += weight; } -#ifdef INFOF - fprintf(stderr,"RM-PF: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step, cr->nodeid, - erg->xc_ref_ind[j],erg->f_rot_loc[j][XX], erg->f_rot_loc[j][YY], erg->f_rot_loc[j][ZZ],erg->fix_torque_v); -#endif + + PRINT_FORCE_J + } /* end of loop over local rotation group atoms */ erg->V = 0.5*rotg->k*V; + + PRINT_POT_TAU } @@ -2628,7 +2618,7 @@ static void do_radial_motion2( rvec x[], /* The positions */ matrix box, /* The simulation box */ double t, /* Time in picoseconds */ - int step, /* The time step */ + gmx_large_int_t step, /* The time step */ gmx_bool bTorque) { int ii,iigrp,j; @@ -2638,7 +2628,7 @@ static void do_radial_motion2( gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ rvec xj_u; /* xj - u */ rvec tmpvec,tmpvec2; - real fac,fac2,Vpart; + real fac,fac2,Vpart=0.0; rvec rj,sj; real sjrj; rvec v_xj_u; /* v x (xj - u) */ @@ -2664,13 +2654,6 @@ static void do_radial_motion2( radial_motion2_precalc_inner_sum(rotg,innersumvec); } - /* Clear values from last time step */ - erg->V = 0.0; - Vpart = 0.0; - erg->fix_torque_v = 0.0; - erg->fix_angles_v = 0.0; - erg->fix_weight_v = 0.0; - N_M = rotg->nat * erg->invmass; /* Each process calculates the forces on its local atoms */ @@ -2743,19 +2726,20 @@ static void do_radial_motion2( if (bTorque) { /* Add to the torque of this rotation group */ - erg->fix_torque_v += torque(rotg->vec, erg->f_rot_loc[j], xj, erg->xc_center); + erg->torque_v += torque(rotg->vec, erg->f_rot_loc[j], xj, erg->xc_center); /* Calculate the angle between reference and actual rotation group atom. */ angle(rotg, xj_u, rj, &alpha, &weight); /* angle in rad, weighted */ - erg->fix_angles_v += alpha * weight; - erg->fix_weight_v += weight; + erg->angle_v += alpha * weight; + erg->weight_v += weight; } -#ifdef INFOF - fprintf(stderr,"RM2: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step, cr->nodeid, - erg->xc_ref_ind[j],erg->f_rot_loc[j][XX], erg->f_rot_loc[j][YY], erg->f_rot_loc[j][ZZ],erg->fix_torque_v); -#endif + + PRINT_FORCE_J + } /* end of loop over local rotation group atoms */ erg->V = 0.5*rotg->k*Vpart; + + PRINT_POT_TAU } @@ -2894,11 +2878,9 @@ static void init_rot_group(FILE *fplog,t_commrec *cr,int g,t_rotgrp *rotg, /* Do we have a flexible axis? */ - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); - + bFlex = ISFLEX(rotg); /* Do we use a global set of coordinates? */ - bColl = bFlex || (rotg->eType==erotgRMPF) || (rotg->eType==erotgRM2PF); + bColl = ISCOLL(rotg); erg=rotg->enfrotgrp; @@ -3072,7 +3054,8 @@ extern void init_rot(FILE *fplog,t_inputrec *ir,int nfile,const t_filenm fnm[], int nat_max=0; /* Size of biggest rotation group */ gmx_enfrot_t er; /* Pointer to the enforced rotation buffer variables */ gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */ - rvec *x_pbc; /* Space for the pbc-correct atom positions */ + rvec *x_pbc=NULL; /* Space for the pbc-correct atom positions */ + gmx_bool bHaveFlexGroups = FALSE; if ( (PAR(cr)) && !DOMAINDECOMP(cr) ) @@ -3092,7 +3075,7 @@ extern void init_rot(FILE *fplog,t_inputrec *ir,int nfile,const t_filenm fnm[], if (fplog) fprintf(fplog, "%s rerun - will write rotation output every available step.\n", RotStr); rot->nstrout = 1; - rot->nsttout = 1; + rot->nstsout = 1; } er->out_slabs = NULL; @@ -3108,9 +3091,11 @@ extern void init_rot(FILE *fplog,t_inputrec *ir,int nfile,const t_filenm fnm[], do_pbc_first_mtop(NULL,ir->ePBC,box,mtop,x_pbc); } - for(g=0; gngrp; g++) + for (g=0; gngrp; g++) { rotg = &rot->grp[g]; + if (ISFLEX(rotg)) + bHaveFlexGroups = TRUE; if (fplog) fprintf(fplog,"%s group %d type '%s'\n", RotStr, g, erotg_names[rotg->eType]); @@ -3156,12 +3141,11 @@ extern void init_rot(FILE *fplog,t_inputrec *ir,int nfile,const t_filenm fnm[], if (MASTER(cr)) { er->out_rot = open_rot_out(opt2fn("-ro",nfile,fnm), rot, oenv, Flags); - if ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ) + if (bHaveFlexGroups) { if (rot->nstrout > 0) er->out_angles = open_angles_out(rot, opt2fn("-ra",nfile,fnm)); - if (rot->nsttout > 0) + if (rot->nstsout > 0) er->out_torque = open_torque_out(rot, opt2fn("-rt",nfile,fnm)); } sfree(x_pbc); @@ -3273,14 +3257,14 @@ extern void do_rotation( matrix box, rvec x[], real t, - int step, + gmx_large_int_t step, gmx_wallcycle_t wcycle, gmx_bool bNS) { int g,i,ii; t_rot *rot; t_rotgrp *rotg; - gmx_bool outstep_torque; + gmx_bool outstep_slab, outstep_rot; gmx_bool bFlex,bColl; float cycles_rot; gmx_enfrot_t er; /* Pointer to the enforced rotation buffer variables */ @@ -3294,11 +3278,13 @@ extern void do_rotation( rot=ir->rot; er=rot->enfrot; - /* At which time steps do we want to output the torque */ - outstep_torque = do_per_step(step, rot->nsttout); - + /* When to output in main rotation output file */ + outstep_rot = do_per_step(step, rot->nstrout); + /* When to output per-slab data */ + outstep_slab = do_per_step(step, rot->nstsout); + /* Output time into rotation output file */ - if (outstep_torque && MASTER(cr)) + if (outstep_rot && MASTER(cr)) fprintf(er->out_rot, "%12.3e",t); /**************************************************************************/ @@ -3309,11 +3295,9 @@ extern void do_rotation( erg=rotg->enfrotgrp; /* Do we have a flexible axis? */ - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); - + bFlex = ISFLEX(rotg); /* Do we use a collective (global) set of coordinates? */ - bColl = bFlex || (rotg->eType==erotgRMPF) || (rotg->eType==erotgRM2PF); + bColl = ISCOLL(rotg); /* Calculate the rotation matrix for this angle: */ erg->degangle = rotg->rate * t; @@ -3369,31 +3353,35 @@ extern void do_rotation( rotg = &rot->grp[g]; erg=rotg->enfrotgrp; - bFlex = ( (rotg->eType==erotgFLEX ) || (rotg->eType==erotgFLEXT ) - || (rotg->eType==erotgFLEX2) || (rotg->eType==erotgFLEX2T) ); + bFlex = ISFLEX(rotg); + bColl = ISCOLL(rotg); - bColl = bFlex || (rotg->eType==erotgRMPF) || (rotg->eType==erotgRM2PF); - - if (outstep_torque && MASTER(cr)) + if (outstep_rot && MASTER(cr)) fprintf(er->out_rot, "%12.4f", erg->degangle); + /* Clear values from last time step */ + erg->V = 0.0; + erg->torque_v = 0.0; + erg->angle_v = 0.0; + erg->weight_v = 0.0; + switch(rotg->eType) { case erotgISO: case erotgISOPF: case erotgPM: case erotgPMPF: - do_fixed(cr,rotg,x,box,t,step,outstep_torque); + do_fixed(cr,rotg,x,box,t,step,outstep_rot); break; case erotgRM: - do_radial_motion(cr,rotg,x,box,t,step,outstep_torque); + do_radial_motion(cr,rotg,x,box,t,step,outstep_rot); break; case erotgRMPF: - do_radial_motion_pf(cr,rotg,x,box,t,step,outstep_torque); + do_radial_motion_pf(cr,rotg,x,box,t,step,outstep_rot); break; case erotgRM2: case erotgRM2PF: - do_radial_motion2(cr,rotg,x,box,t,step,outstep_torque); + do_radial_motion2(cr,rotg,x,box,t,step,outstep_rot); break; case erotgFLEXT: case erotgFLEX2T: @@ -3403,13 +3391,13 @@ extern void do_rotation( get_center(erg->xc, erg->mc, rotg->nat, erg->xc_center); svmul(-1.0, erg->xc_center, transvec); translate_x(erg->xc, rotg->nat, transvec); - do_flexible(cr,er,rotg,g,x,box,t,step,outstep_torque); + do_flexible(cr,er,rotg,g,x,box,t,step,outstep_rot,outstep_slab); break; case erotgFLEX: case erotgFLEX2: /* Do NOT subtract the center of mass in the low level routines! */ clear_rvec(erg->xc_center); - do_flexible(cr,er,rotg,g,x,box,t,step,outstep_torque); + do_flexible(cr,er,rotg,g,x,box,t,step,outstep_rot,outstep_slab); break; default: gmx_fatal(FARGS, "No such rotation potential."); diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c index b176ccf6af..5184758e9a 100644 --- a/src/mdlib/sim_util.c +++ b/src/mdlib/sim_util.c @@ -854,7 +854,7 @@ void do_force(FILE *fplog,t_commrec *cr, /* Add the forces from enforced rotation potentials (if any) */ if (inputrec->bRot) - enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr,step,t); + enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr, step, t); if (PAR(cr) && !(cr->duty & DUTY_PME)) -- 2.22.0