From: Roland Schulz Date: Tue, 26 Mar 2013 19:15:20 +0000 (-0400) Subject: Merge release-4-5-patches into release-4-6 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=d4e96c5f971f37566e159e3b4b32cb189de77714;p=alexxy%2Fgromacs.git Merge release-4-5-patches into release-4-6 Conflicts: CMakeLists.txt configure.ac src/tools/gmx_mindist.c Change-Id: I303318f17527a8e1bfe494eb174565b8d5491870 --- d4e96c5f971f37566e159e3b4b32cb189de77714 diff --cc src/tools/gmx_mindist.c index f86f22ac99,dc9ae1f74f..c44c181412 --- a/src/tools/gmx_mindist.c +++ b/src/tools/gmx_mindist.c @@@ -60,515 -57,386 +60,515 @@@ #include "gmx_ana.h" -static void periodic_dist(matrix box,rvec x[],int n,atom_id index[], - real *rmin,real *rmax,int *min_ind) +static void periodic_dist(matrix box, rvec x[], int n, atom_id index[], + real *rmin, real *rmax, int *min_ind) { #define NSHIFT 26 - int sx,sy,sz,i,j,s; - real sqr_box,r2min,r2max,r2; - rvec shift[NSHIFT],d0,d; - - sqr_box = sqr(min(norm(box[XX]),min(norm(box[YY]),norm(box[ZZ])))); - - s = 0; - for(sz=-1; sz<=1; sz++) - for(sy=-1; sy<=1; sy++) - for(sx=-1; sx<=1; sx++) - if (sx!=0 || sy!=0 || sz!=0) { - for(i=0; i r2max) - r2max = r2; - for(s=0; s r2max) + { + r2max = r2; + } + for (s = 0; s < NSHIFT; s++) + { + rvec_add(d0, shift[s], d); + r2 = norm2(d); + if (r2 < r2min) + { + r2min = r2; + min_ind[0] = i; + min_ind[1] = j; + } + } + } + } + + *rmin = sqrt(r2min); + *rmax = sqrt(r2max); } -static void periodic_mindist_plot(const char *trxfn,const char *outfn, - t_topology *top,int ePBC, - int n,atom_id index[],gmx_bool bSplit, +static void periodic_mindist_plot(const char *trxfn, const char *outfn, + t_topology *top, int ePBC, + int n, atom_id index[], gmx_bool bSplit, const output_env_t oenv) { - FILE *out; - const char *leg[5] = { "min per.","max int.","box1","box2","box3" }; - t_trxstatus *status; - real t; - rvec *x; - matrix box; - int natoms,ind_min[2]={0,0},ind_mini=0,ind_minj=0; - real r,rmin,rmax,rmint,tmint; - gmx_bool bFirst; - gmx_rmpbc_t gpbc=NULL; - - natoms=read_first_x(oenv,&status,trxfn,&t,&x,box); - - check_index(NULL,n,index,NULL,natoms); - - out = xvgropen(outfn,"Minimum distance to periodic image", - output_env_get_time_label(oenv),"Distance (nm)",oenv); - if (output_env_get_print_xvgr_codes(oenv)) - fprintf(out,"@ subtitle \"and maximum internal distance\"\n"); - xvgr_legend(out,5,leg,oenv); - - rmint = box[XX][XX]; - tmint = 0; - - if (NULL != top) - gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); - - bFirst=TRUE; - do { - if (NULL != top) - gmx_rmpbc(gpbc,natoms,box,x); - - periodic_dist(box,x,n,index,&rmin,&rmax,ind_min); - if (rmin < rmint) { - rmint = rmin; - tmint = t; - ind_mini = ind_min[0]; - ind_minj = ind_min[1]; - } - if ( bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv))<1e-5 ) - fprintf(out, "&\n"); - fprintf(out,"\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n", - output_env_conv_time(oenv,t),rmin,rmax,norm(box[0]),norm(box[1]),norm(box[2])); - bFirst=FALSE; - } while(read_next_x(oenv,status,&t,natoms,x,box)); - - if (NULL != top) - gmx_rmpbc_done(gpbc); - - ffclose(out); - - fprintf(stdout, - "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n" - "between atoms %d and %d\n", - rmint,output_env_conv_time(oenv,tmint),output_env_get_time_unit(oenv), - index[ind_mini]+1,index[ind_minj]+1); + FILE *out; + const char *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" }; + t_trxstatus *status; + real t; + rvec *x; + matrix box; + int natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0; + real r, rmin, rmax, rmint, tmint; + gmx_bool bFirst; + gmx_rmpbc_t gpbc = NULL; + + natoms = read_first_x(oenv, &status, trxfn, &t, &x, box); + + check_index(NULL, n, index, NULL, natoms); + + out = xvgropen(outfn, "Minimum distance to periodic image", + output_env_get_time_label(oenv), "Distance (nm)", oenv); + if (output_env_get_print_xvgr_codes(oenv)) + { + fprintf(out, "@ subtitle \"and maximum internal distance\"\n"); + } + xvgr_legend(out, 5, leg, oenv); + + rmint = box[XX][XX]; + tmint = 0; + + if (NULL != top) + { + gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box); + } + + bFirst = TRUE; + do + { + if (NULL != top) + { + gmx_rmpbc(gpbc, natoms, box, x); + } + + periodic_dist(box, x, n, index, &rmin, &rmax, ind_min); + if (rmin < rmint) + { + rmint = rmin; + tmint = t; + ind_mini = ind_min[0]; + ind_minj = ind_min[1]; + } + if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5) + { + fprintf(out, "&\n"); + } + fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n", + output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2])); + bFirst = FALSE; + } + while (read_next_x(oenv, status, &t, natoms, x, box)); + + if (NULL != top) + { + gmx_rmpbc_done(gpbc); + } + + ffclose(out); + + fprintf(stdout, + "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n" + "between atoms %d and %d\n", + rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv), + index[ind_mini]+1, index[ind_minj]+1); } -static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[], - int nx1,int nx2, atom_id index1[], atom_id index2[], - gmx_bool bGroup, - real *rmin, real *rmax, int *nmin, int *nmax, - int *ixmin, int *jxmin, int *ixmax, int *jxmax) +static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[], + int nx1, int nx2, atom_id index1[], atom_id index2[], + gmx_bool bGroup, + real *rmin, real *rmax, int *nmin, int *nmax, + int *ixmin, int *jxmin, int *ixmax, int *jxmax) { - int i,j,i0=0,j1; - int ix,jx; - atom_id *index3; - rvec dx; - real r2,rmin2,rmax2,rcut2; - t_pbc pbc; - int nmin_j,nmax_j; - - *ixmin = -1; - *jxmin = -1; - *ixmax = -1; - *jxmax = -1; - *nmin = 0; - *nmax = 0; - - rcut2=sqr(rcut); - - /* Must init pbc every step because of pressure coupling */ - if (bPBC) - set_pbc(&pbc,ePBC,box); - if (index2) { - i0=0; - j1=nx2; - index3=index2; - } else { - j1=nx1; - index3=index1; - } - - rmin2=1e12; - rmax2=-1e12; - - for(j=0; (j < j1); j++) { - jx = index3[j]; - if (index2 == NULL) { - i0 = j + 1; - } - nmin_j = 0; - nmax_j = 0; - for(i=i0; (i < nx1); i++) { - ix = index1[i]; - if (ix != jx) { - if (bPBC) - pbc_dx(&pbc,x[ix],x[jx],dx); - else - rvec_sub(x[ix],x[jx],dx); - r2=iprod(dx,dx); - if (r2 < rmin2) { - rmin2=r2; - *ixmin=ix; - *jxmin=jx; - } - if (r2 > rmax2) { - rmax2=r2; - *ixmax=ix; - *jxmax=jx; - } - if (r2 <= rcut2) { - nmin_j++; - } else if (r2 > rcut2) { - nmax_j++; - } - } - } - if (bGroup) { - if (nmin_j > 0) { - (*nmin)++; - } - if (nmax_j > 0) { - (*nmax)++; - } - } else { - *nmin += nmin_j; - *nmax += nmax_j; - } - } - *rmin = sqrt(rmin2); - *rmax = sqrt(rmax2); + int i, j, i0 = 0, j1; + int ix, jx; + atom_id *index3; + rvec dx; + real r2, rmin2, rmax2, rcut2; + t_pbc pbc; + int nmin_j, nmax_j; + + *ixmin = -1; + *jxmin = -1; + *ixmax = -1; + *jxmax = -1; + *nmin = 0; + *nmax = 0; + + rcut2 = sqr(rcut); + + /* Must init pbc every step because of pressure coupling */ + if (bPBC) + { + set_pbc(&pbc, ePBC, box); + } + if (index2) + { + i0 = 0; + j1 = nx2; + index3 = index2; + } + else + { + j1 = nx1; + index3 = index1; + } + + rmin2 = 1e12; + rmax2 = -1e12; + + for (j = 0; (j < j1); j++) + { + jx = index3[j]; + if (index2 == NULL) + { + i0 = j + 1; + } + nmin_j = 0; + nmax_j = 0; + for (i = i0; (i < nx1); i++) + { + ix = index1[i]; + if (ix != jx) + { + if (bPBC) + { + pbc_dx(&pbc, x[ix], x[jx], dx); + } + else + { + rvec_sub(x[ix], x[jx], dx); + } + r2 = iprod(dx, dx); + if (r2 < rmin2) + { + rmin2 = r2; + *ixmin = ix; + *jxmin = jx; + } + if (r2 > rmax2) + { + rmax2 = r2; + *ixmax = ix; + *jxmax = jx; + } + if (r2 <= rcut2) + { + nmin_j++; + } + else if (r2 > rcut2) + { + nmax_j++; + } + } + } + if (bGroup) + { + if (nmin_j > 0) + { + (*nmin)++; + } + if (nmax_j > 0) + { + (*nmax)++; + } + } + else + { + *nmin += nmin_j; + *nmax += nmax_j; + } + } + *rmin = sqrt(rmin2); + *rmax = sqrt(rmax2); } -void dist_plot(const char *fn,const char *afile,const char *dfile, - const char *nfile,const char *rfile,const char *xfile, - real rcut,gmx_bool bMat,t_atoms *atoms, - int ng,atom_id *index[],int gnx[],char *grpn[],gmx_bool bSplit, - gmx_bool bMin, int nres, atom_id *residue,gmx_bool bPBC,int ePBC, - gmx_bool bGroup,gmx_bool bEachResEachTime, gmx_bool bPrintResName, +void dist_plot(const char *fn, const char *afile, const char *dfile, + const char *nfile, const char *rfile, const char *xfile, + real rcut, gmx_bool bMat, t_atoms *atoms, + int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit, + gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC, + gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName, const output_env_t oenv) { - FILE *atm,*dist,*num; - t_trxstatus *trxout; - char buf[256]; - char **leg; - real t,dmin,dmax,**mindres=NULL,**maxdres=NULL; - int nmin,nmax; - t_trxstatus *status; - int i=-1,j,k,natoms; - int min1,min2,max1,max2,min1r,min2r,max1r,max2r; - atom_id oindex[2]; - rvec *x0; - matrix box; - t_trxframe frout; - gmx_bool bFirst; - FILE *respertime=NULL; - - if ((natoms=read_first_x(oenv,&status,fn,&t,&x0,box))==0) - gmx_fatal(FARGS,"Could not read coordinates from statusfile\n"); - - sprintf(buf,"%simum Distance",bMin ? "Min" : "Max"); - dist= xvgropen(dfile,buf,output_env_get_time_label(oenv),"Distance (nm)",oenv); - sprintf(buf,"Number of Contacts %s %g nm",bMin ? "<" : ">",rcut); - num = nfile ? xvgropen(nfile,buf,output_env_get_time_label(oenv),"Number",oenv) : NULL; - atm = afile ? ffopen(afile,"w") : NULL; - trxout = xfile ? open_trx(xfile,"w") : NULL; - - if (bMat) { - if (ng == 1) { - snew(leg,1); - sprintf(buf,"Internal in %s",grpn[0]); - leg[0]=strdup(buf); - xvgr_legend(dist,0,(const char**)leg,oenv); - if (num) xvgr_legend(num,0,(const char**)leg,oenv); - } - else { - snew(leg,(ng*(ng-1))/2); - for(i=j=0; (iresinfo[atoms->atom[index[0][residue[j]]].resind].name),atoms->atom[index[0][residue[j]]].resind); - fprintf(respertime,"\n"); - - } - - j=0; - if (nres) { - snew(mindres, ng-1); - snew(maxdres, ng-1); - for(i=1; i", rcut); + num = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL; + atm = afile ? ffopen(afile, "w") : NULL; + trxout = xfile ? open_trx(xfile, "w") : NULL; + + if (bMat) + { + if (ng == 1) + { + snew(leg, 1); + sprintf(buf, "Internal in %s", grpn[0]); + leg[0] = strdup(buf); + xvgr_legend(dist, 0, (const char**)leg, oenv); + if (num) + { + xvgr_legend(num, 0, (const char**)leg, oenv); + } + } + else + { + snew(leg, (ng*(ng-1))/2); + for (i = j = 0; (i < ng-1); i++) + { + for (k = i+1; (k < ng); k++, j++) + { + sprintf(buf, "%s-%s", grpn[i], grpn[k]); + leg[j] = strdup(buf); + } + } + xvgr_legend(dist, j, (const char**)leg, oenv); + if (num) + { + xvgr_legend(num, j, (const char**)leg, oenv); + } + } + } + else + { + snew(leg, ng-1); + for (i = 0; (i < ng-1); i++) + { + sprintf(buf, "%s-%s", grpn[0], grpn[i+1]); + leg[i] = strdup(buf); + } + xvgr_legend(dist, ng-1, (const char**)leg, oenv); + if (num) + { + xvgr_legend(num, ng-1, (const char**)leg, oenv); + } + } + + if (bEachResEachTime) + { + sprintf(buf, "%simum Distance", bMin ? "Min" : "Max"); + respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv); + xvgr_legend(respertime, ng-1, (const char**)leg, oenv); + if (bPrintResName) + { + fprintf(respertime, "# "); + } + for (j = 0; j < nres; j++) + { + fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind); + } + fprintf(respertime, "\n"); + + } + + j = 0; + if (nres) + { + snew(mindres, ng-1); + snew(maxdres, ng-1); + for (i = 1; i < ng; i++) + { + snew(mindres[i-1], nres); + snew(maxdres[i-1], nres); + for (j = 0; j < nres; j++) + { + mindres[i-1][j] = 1e6; + } + /* maxdres[*][*] is already 0 */ + } + } + bFirst = TRUE; + do + { + if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5) + { + fprintf(dist, "&\n"); + if (num) + { + fprintf(num, "&\n"); + } + if (atm) + { + fprintf(atm, "&\n"); + } + } + fprintf(dist, "%12e", output_env_conv_time(oenv, t)); + if (num) + { + fprintf(num, "%12e", output_env_conv_time(oenv, t)); + } + + if (bMat) + { + if (ng == 1) + { + calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup, + &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2); + fprintf(dist, " %12e", bMin ? dmin : dmax); + if (num) + { + fprintf(num, " %8d", bMin ? nmin : nmax); + } + } + else + { + for (i = 0; (i < ng-1); i++) + { + for (k = i+1; (k < ng); k++) + { + calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k], + bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2); + fprintf(dist, " %12e", bMin ? dmin : dmax); + if (num) + { + fprintf(num, " %8d", bMin ? nmin : nmax); + } + } + } + } + } + else + { + for (i = 1; (i < ng); i++) + { + calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup, + &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2); + fprintf(dist, " %12e", bMin ? dmin : dmax); + if (num) + { + fprintf(num, " %8d", bMin ? nmin : nmax); + } + if (nres) + { + for (j = 0; j < nres; j++) + { + calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i], + &(index[0][residue[j]]), index[i], bGroup, + &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r); + mindres[i-1][j] = min(mindres[i-1][j], dmin); + maxdres[i-1][j] = max(maxdres[i-1][j], dmax); + } + } + } + } + fprintf(dist, "\n"); + if (num) + { + fprintf(num, "\n"); + } + if ( (bMin ? min1 : max1) != -1) + { + if (atm) + { + fprintf(atm, "%12e %12d %12d\n", + output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1), + 1+(bMin ? min2 : max2)); + } + } + + if (trxout) + { + oindex[0] = bMin ? min1 : max1; + oindex[1] = bMin ? min2 : max2; + write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL); + } + bFirst = FALSE; + /*dmin should be minimum distance for residue and group*/ + if (bEachResEachTime) + { + fprintf(respertime, "%12e", t); + for (i = 1; i < ng; i++) + { + for (j = 0; j < nres; j++) + { + fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]); + /*reset distances for next time point*/ + mindres[i-1][j] = 1e6; + maxdres[i-1][j] = 0; + } + } + fprintf(respertime, "\n"); + } + } + while (read_next_x(oenv, status, &t, natoms, x0, box)); + + close_trj(status); + ffclose(dist); + if (num) + { + ffclose(num); + } + if (atm) + { + ffclose(atm); + } + if (trxout) + { + close_trx(trxout); + } + + if (nres && !bEachResEachTime) + { + FILE *res; + + sprintf(buf, "%simum Distance", bMin ? "Min" : "Max"); + res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv); + xvgr_legend(res, ng-1, (const char**)leg, oenv); + for (j = 0; j < nres; j++) + { + fprintf(res, "%4d", j+1); + for (i = 1; i < ng; i++) + { + fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]); + } + fprintf(res, "\n"); + } + } + + sfree(x0); } int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)