From aef61d9804d107c75e373f3e4c29db5c8f18d6a4 Mon Sep 17 00:00:00 2001 From: Erik Lindahl Date: Fri, 10 Jul 2015 01:27:57 +0200 Subject: [PATCH] Convert support files in gmxana to C++ Renamed symbols using reserved C++ names, and removed geminate and ballistic code paths that just resulted in GROMACS errors messages about missing code when they were called. Unused parts of gmx_hbond have been cleaned up according to comments from Erik Marklund. Change-Id: I925b1b5ca4a8420e0c8a38ca99f56361359d713b --- src/gromacs/gmxana/{anadih.c => anadih.cpp} | 62 +- .../gmxana/{binsearch.c => binsearch.cpp} | 38 +- src/gromacs/gmxana/binsearch.h | 16 +- src/gromacs/gmxana/{cmat.c => cmat.cpp} | 8 +- src/gromacs/gmxana/cmat.h | 11 +- .../gmxana/{dens_filter.c => dens_filter.cpp} | 29 +- src/gromacs/gmxana/{dlist.c => dlist.cpp} | 87 +- src/gromacs/gmxana/{edittop.c => edittop.cpp} | 6 +- src/gromacs/gmxana/{eigio.c => eigio.cpp} | 0 src/gromacs/gmxana/eigio.h | 9 + src/gromacs/gmxana/{fitahx.c => fitahx.cpp} | 12 +- src/gromacs/gmxana/fitahx.h | 15 +- src/gromacs/gmxana/geminate.c | 767 ------- src/gromacs/gmxana/geminate.h | 184 -- src/gromacs/gmxana/gmx_analyze.c | 163 +- src/gromacs/gmxana/gmx_hbond.c | 1920 +++-------------- src/gromacs/gmxana/{hxprops.c => hxprops.cpp} | 41 +- src/gromacs/gmxana/hxprops.h | 11 +- src/gromacs/gmxana/{nrama.c => nrama.cpp} | 4 +- .../gmxana/{nsfactor.c => nsfactor.cpp} | 36 +- .../gmxana/{powerspect.c => powerspect.cpp} | 2 +- .../gmxana/{pp2shift.c => pp2shift.cpp} | 12 +- src/gromacs/gmxana/{princ.c => princ.cpp} | 6 +- src/gromacs/gmxana/{sfactor.c => sfactor.cpp} | 48 +- 24 files changed, 541 insertions(+), 2946 deletions(-) rename src/gromacs/gmxana/{anadih.c => anadih.cpp} (94%) rename src/gromacs/gmxana/{binsearch.c => binsearch.cpp} (88%) rename src/gromacs/gmxana/{cmat.c => cmat.cpp} (97%) rename src/gromacs/gmxana/{dens_filter.c => dens_filter.cpp} (88%) rename src/gromacs/gmxana/{dlist.c => dlist.cpp} (81%) rename src/gromacs/gmxana/{edittop.c => edittop.cpp} (98%) rename src/gromacs/gmxana/{eigio.c => eigio.cpp} (100%) rename src/gromacs/gmxana/{fitahx.c => fitahx.cpp} (94%) delete mode 100644 src/gromacs/gmxana/geminate.c delete mode 100644 src/gromacs/gmxana/geminate.h rename src/gromacs/gmxana/{hxprops.c => hxprops.cpp} (94%) rename src/gromacs/gmxana/{nrama.c => nrama.cpp} (99%) rename src/gromacs/gmxana/{nsfactor.c => nsfactor.cpp} (86%) rename src/gromacs/gmxana/{powerspect.c => powerspect.cpp} (98%) rename src/gromacs/gmxana/{pp2shift.c => pp2shift.cpp} (96%) rename src/gromacs/gmxana/{princ.c => princ.cpp} (98%) rename src/gromacs/gmxana/{sfactor.c => sfactor.cpp} (93%) diff --git a/src/gromacs/gmxana/anadih.c b/src/gromacs/gmxana/anadih.cpp similarity index 94% rename from src/gromacs/gmxana/anadih.c rename to src/gromacs/gmxana/anadih.cpp index a1fc18e8d1..43b50af15e 100644 --- a/src/gromacs/gmxana/anadih.c +++ b/src/gromacs/gmxana/anadih.cpp @@ -36,9 +36,11 @@ */ #include "gmxpre.h" -#include -#include -#include +#include +#include +#include + +#include #include "gromacs/fileio/confio.h" #include "gromacs/fileio/trxio.h" @@ -167,7 +169,7 @@ void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans, char title[256]; int i, j, k, Dih, ntrans; int cur_bin, new_bin; - real ttime, tt; + real ttime; real *rot_occ[NROT]; int (*calc_bin)(real, int, real); real dt; @@ -243,8 +245,8 @@ void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans, prev = dih[i][j]; - mind = min(mind, dih[i][j]); - maxd = max(maxd, dih[i][j]); + mind = std::min(mind, dih[i][j]); + maxd = std::max(maxd, dih[i][j]); if ( (maxd - mind) > 2*M_PI/3) /* or 120 degrees, assuming */ { /* multiplicity 3. Not so general.*/ tr_f[j]++; @@ -357,7 +359,7 @@ void mk_multiplicity_lookup (int *multiplicity, int maxchi, { for (i = 0; (i < nlist); i++) { - strncpy(name, dlist[i].name, 3); + std::strncpy(name, dlist[i].name, 3); name[3] = '\0'; if (((Dih < edOmega) ) || ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) || @@ -375,16 +377,16 @@ void mk_multiplicity_lookup (int *multiplicity, int maxchi, /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/ if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1)) { - if ( ((strstr(name, "PHE") != NULL) && (Dih == edChi2)) || - ((strstr(name, "TYR") != NULL) && (Dih == edChi2)) || - ((strstr(name, "PTR") != NULL) && (Dih == edChi2)) || - ((strstr(name, "TRP") != NULL) && (Dih == edChi2)) || - ((strstr(name, "HIS") != NULL) && (Dih == edChi2)) || - ((strstr(name, "GLU") != NULL) && (Dih == edChi3)) || - ((strstr(name, "ASP") != NULL) && (Dih == edChi2)) || - ((strstr(name, "GLN") != NULL) && (Dih == edChi3)) || - ((strstr(name, "ASN") != NULL) && (Dih == edChi2)) || - ((strstr(name, "ARG") != NULL) && (Dih == edChi4)) ) + if ( ((std::strstr(name, "PHE") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "TYR") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "PTR") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "TRP") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "HIS") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "GLU") != NULL) && (Dih == edChi3)) || + ((std::strstr(name, "ASP") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "GLN") != NULL) && (Dih == edChi3)) || + ((std::strstr(name, "ASN") != NULL) && (Dih == edChi2)) || + ((std::strstr(name, "ARG") != NULL) && (Dih == edChi4)) ) { multiplicity[j] = 2; } @@ -510,7 +512,6 @@ void get_chi_product_traj (real **dih, int nframes, int nlist, index = lookup[i][0]; /* index into dih of chi1 of res i */ if (index == -1) { - b = 0; bRotZero = TRUE; bHaveChi = FALSE; } @@ -652,15 +653,14 @@ void calc_distribution_props(int nh, int histo[], real start, { d = invth*histo[j]; ang = j*fac-start; - c1 = cos(ang); - c2 = c1*c1; + c1 = std::cos(ang); dc = d*c1; - ds = d*sin(ang); + ds = d*std::sin(ang); tdc += dc; tds += ds; for (i = 0; (i < nkkk); i++) { - c1 = cos(ang+kkk[i].offset); + c1 = std::cos(ang+kkk[i].offset); c2 = c1*c1; Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C); kkk[i].Jc += histo[j]*Jc; @@ -670,7 +670,7 @@ void calc_distribution_props(int nh, int histo[], real start, for (i = 0; (i < nkkk); i++) { kkk[i].Jc /= th; - kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc)); + kkk[i].Jcsig = std::sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc)); } *S2 = tdc*tdc+tds*tds; } @@ -680,7 +680,7 @@ static void calc_angles(struct t_pbc *pbc, { int i, ix, t1, t2; rvec r_ij, r_kj; - real costh; + real costh = 0.0; for (i = ix = 0; (ix < n3); i++, ix += 3) { @@ -759,8 +759,8 @@ void make_histo(FILE *log, minx = maxx = data[0]; for (i = 1; (i < ndata); i++) { - minx = min(minx, data[i]); - maxx = max(maxx, data[i]); + minx = std::min(minx, data[i]); + maxx = std::max(maxx, data[i]); } fprintf(log, "Min data: %10g Max data: %10g\n", minx, maxx); } @@ -773,7 +773,7 @@ void make_histo(FILE *log, } for (i = 0; (i < ndata); i++) { - ind = (data[i]-minx)*dx; + ind = static_cast((data[i]-minx)*dx); if ((ind >= 0) && (ind < npoints)) { histo[ind]++; @@ -819,7 +819,7 @@ void read_ang_dih(const char *trj_fn, { struct t_pbc *pbc; t_trxstatus *status; - int i, angind, natoms, total, teller; + int i, angind, total, teller; int nangles, n_alloc; real t, fraction, pifac, aa, angle; real *angles[2]; @@ -829,7 +829,7 @@ void read_ang_dih(const char *trj_fn, #define prev (1-cur) snew(pbc, 1); - natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box); + read_first_x(oenv, &status, trj_fn, &t, &x, box); if (bAngles) { @@ -911,7 +911,7 @@ void read_ang_dih(const char *trj_fn, for (i = 0; (i < nangles); i++) { real dd = angles[cur][i]; - angles[cur][i] = atan2(sin(dd), cos(dd)); + angles[cur][i] = std::atan2(sin(dd), cos(dd)); } } else @@ -961,7 +961,7 @@ void read_ang_dih(const char *trj_fn, } /* Update the distribution histogram */ - angind = (int) ((angle*maxangstat)/pifac + 0.5); + angind = static_cast((angle*maxangstat)/pifac + 0.5); if (angind == maxangstat) { angind = 0; diff --git a/src/gromacs/gmxana/binsearch.c b/src/gromacs/gmxana/binsearch.cpp similarity index 88% rename from src/gromacs/gmxana/binsearch.c rename to src/gromacs/gmxana/binsearch.cpp index ac5135bef4..1c9f474d14 100644 --- a/src/gromacs/gmxana/binsearch.c +++ b/src/gromacs/gmxana/binsearch.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -34,7 +34,9 @@ */ #include "gmxpre.h" -#include +#include "binsearch.h" + +#include #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/real.h" @@ -109,44 +111,44 @@ void insertionSort(real *arr, int *perm, int startndx, int endndx, int direction int BinarySearch (real *array, int low, int high, real key, int direction) { - int mid, max, min; - max = high+2; - min = low+1; + int iMid, iMax, iMin; + iMax = high+2; + iMin = low+1; /*Iterative implementation*/ if (direction >= 0) { - while (max-min > 1) + while (iMax-iMin > 1) { - mid = (min+max)>>1; - if (key < array[mid-1]) + iMid = (iMin+iMax)>>1; + if (key < array[iMid-1]) { - max = mid; + iMax = iMid; } else { - min = mid; + iMin = iMid; } } - return min; + return iMin; } else if (direction < 0) { - while (max-min > 1) + while (iMax-iMin > 1) { - mid = (min+max)>>1; - if (key > array[mid-1]) + iMid = (iMin+iMax)>>1; + if (key > array[iMid-1]) { - max = mid; + iMax = iMid; } else { - min = mid; + iMin = iMid; } } - return min-1; + return iMin-1; } /*end -ifelse direction*/ return -1; @@ -169,7 +171,6 @@ int LinearSearch (double *array, int startindx, int stopindx, if (direction >= 0) { - keyindex = startindx; for (i = startindx; i <= stopindx; i++) { (*count)++; @@ -182,7 +183,6 @@ int LinearSearch (double *array, int startindx, int stopindx, } else if (direction < 0) { - keyindex = stopindx; for (i = stopindx; i >= startindx; i--) { (*count)++; diff --git a/src/gromacs/gmxana/binsearch.h b/src/gromacs/gmxana/binsearch.h index 85bb9ec93c..24cf8f465b 100644 --- a/src/gromacs/gmxana/binsearch.h +++ b/src/gromacs/gmxana/binsearch.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -42,17 +42,17 @@ extern "C" { #endif -extern void rangeArray(int *ar, int size); +void rangeArray(int *ar, int size); -extern void insertionSort(real *ar, int *perm, int start, int end, int direction); +void insertionSort(real *ar, int *perm, int start, int end, int direction); -extern int BinarySearch(real *ar, int start, int end, real key, int direction); +int BinarySearch(real *ar, int start, int end, real key, int direction); -extern int start_binsearch(real *array, int *perm, int low, int high, - real key, int direction); +int start_binsearch(real *array, int *perm, int low, int high, + real key, int direction); -extern int LinearSearch(double *array, int startindx, int stopindx, - double key, int *count, int direction); +int LinearSearch(double *array, int startindx, int stopindx, + double key, int *count, int direction); #ifdef __cplusplus } diff --git a/src/gromacs/gmxana/cmat.c b/src/gromacs/gmxana/cmat.cpp similarity index 97% rename from src/gromacs/gmxana/cmat.c rename to src/gromacs/gmxana/cmat.cpp index fab0b3444a..ff9aa87907 100644 --- a/src/gromacs/gmxana/cmat.c +++ b/src/gromacs/gmxana/cmat.cpp @@ -38,6 +38,8 @@ #include "cmat.h" +#include + #include "gromacs/fileio/matio.h" #include "gromacs/fileio/xvgr.h" #include "gromacs/legacyheaders/macros.h" @@ -127,13 +129,13 @@ void reset_index(t_mat *m) void set_mat_entry(t_mat *m, int i, int j, real val) { m->mat[i][j] = m->mat[j][i] = val; - m->maxrms = max(m->maxrms, val); + m->maxrms = std::max(m->maxrms, val); if (j != i) { - m->minrms = min(m->minrms, val); + m->minrms = std::min(m->minrms, val); } m->sumrms += val; - m->nn = max(m->nn, max(j+1, i+1)); + m->nn = std::max(m->nn, std::max(j+1, i+1)); } void done_mat(t_mat **m) diff --git a/src/gromacs/gmxana/cmat.h b/src/gromacs/gmxana/cmat.h index 64673ec135..63f6ca932d 100644 --- a/src/gromacs/gmxana/cmat.h +++ b/src/gromacs/gmxana/cmat.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -40,6 +40,11 @@ #include "gromacs/legacyheaders/typedefs.h" +#ifdef __cplusplus +extern "C" +{ +#endif + typedef struct { int i, j; real dist; @@ -86,4 +91,8 @@ extern void rmsd_distribution(const char *fn, t_mat *m, const output_env_t oenv) extern t_clustid *new_clustid(int n1); +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/gromacs/gmxana/dens_filter.c b/src/gromacs/gmxana/dens_filter.cpp similarity index 88% rename from src/gromacs/gmxana/dens_filter.c rename to src/gromacs/gmxana/dens_filter.cpp index bd08b5247a..3e2f6b2602 100644 --- a/src/gromacs/gmxana/dens_filter.c +++ b/src/gromacs/gmxana/dens_filter.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2011,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2011,2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -41,29 +41,12 @@ #include "dens_filter.h" -#include +#include #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/math/vec.h" #include "gromacs/utility/smalloc.h" -#ifdef GMX_DOUBLE -#define EXP(x) (exp(x)) -#else -#define EXP(x) (expf(x)) -#endif - -/* Modulo function assuming y>0 but with arbitrary INTEGER x */ -static int MOD(int x, int y) -{ - if (x < 0) - { - x = x+y; - } - return (mod(x, y)); -} - - gmx_bool convolution(int dataSize, real *x, int kernelSize, real* kernel) { int i, j, k; @@ -110,7 +93,7 @@ gmx_bool convolution(int dataSize, real *x, int kernelSize, real* kernel) gmx_bool periodic_convolution(int datasize, real *x, int kernelsize, real *kernel) { - int i, j, k, nj; + int i, j, idx; real *filtered; if (!x || !kernel) @@ -128,7 +111,9 @@ gmx_bool periodic_convolution(int datasize, real *x, int kernelsize, { for (j = 0; (j < kernelsize); j++) { - filtered[i] += kernel[j]*x[MOD((i-j), datasize)]; + // add datasize in case i-j is <0 + idx = i-j + datasize; + filtered[i] += kernel[j]*x[idx % datasize]; } } for (i = 0; i < datasize; i++) @@ -153,7 +138,7 @@ void gausskernel(real *out, int n, real var) for (i = -k; i <= k; i++) { arg = (i*i)/(2*var); - tot += out[j++] = EXP(-arg); + tot += out[j++] = std::exp(-arg); } for (i = 0; i < j; i++) { diff --git a/src/gromacs/gmxana/dlist.c b/src/gromacs/gmxana/dlist.cpp similarity index 81% rename from src/gromacs/gmxana/dlist.c rename to src/gromacs/gmxana/dlist.cpp index b7578cef62..1d4f7aaf64 100644 --- a/src/gromacs/gmxana/dlist.c +++ b/src/gromacs/gmxana/dlist.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -36,8 +36,8 @@ */ #include "gmxpre.h" -#include -#include +#include +#include #include "gromacs/gmxana/gstat.h" #include "gromacs/topology/residuetypes.h" @@ -49,7 +49,7 @@ t_dlist *mk_dlist(FILE *log, gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi, int maxchi, int r0, gmx_residuetype_t *rt) { - int ires, i, j, k, ii; + int i, j, ii; t_dihatms atm, prev; int nl = 0, nc[edMax]; char *thisres; @@ -61,11 +61,10 @@ t_dlist *mk_dlist(FILE *log, { nc[i] = 0; } - ires = -1; i = 0; while (i < atoms->nr) { - ires = atoms->atom[i].resind; + int ires = atoms->atom[i].resind; /* Initiate all atom numbers to -1 */ atm.minC = atm.H = atm.N = atm.C = atm.O = atm.minCalpha = -1; @@ -78,71 +77,71 @@ t_dlist *mk_dlist(FILE *log, /* maybe should allow for chis to hydrogens? */ while ((i < atoms->nr) && (atoms->atom[i].resind == ires)) { - if ((strcmp(*(atoms->atomname[i]), "H") == 0) || - (strcmp(*(atoms->atomname[i]), "H1") == 0) || - (strcmp(*(atoms->atomname[i]), "HN") == 0) ) + if ((std::strcmp(*(atoms->atomname[i]), "H") == 0) || + (std::strcmp(*(atoms->atomname[i]), "H1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "HN") == 0) ) { atm.H = i; } - else if (strcmp(*(atoms->atomname[i]), "N") == 0) + else if (std::strcmp(*(atoms->atomname[i]), "N") == 0) { atm.N = i; } - else if (strcmp(*(atoms->atomname[i]), "C") == 0) + else if (std::strcmp(*(atoms->atomname[i]), "C") == 0) { atm.C = i; } - else if ((strcmp(*(atoms->atomname[i]), "O") == 0) || - (strcmp(*(atoms->atomname[i]), "O1") == 0) || - (strcmp(*(atoms->atomname[i]), "OC1") == 0) || - (strcmp(*(atoms->atomname[i]), "OT1") == 0)) + else if ((std::strcmp(*(atoms->atomname[i]), "O") == 0) || + (std::strcmp(*(atoms->atomname[i]), "O1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OC1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OT1") == 0)) { atm.O = i; } - else if (strcmp(*(atoms->atomname[i]), "CA") == 0) + else if (std::strcmp(*(atoms->atomname[i]), "CA") == 0) { atm.Cn[1] = i; } - else if (strcmp(*(atoms->atomname[i]), "CB") == 0) + else if (std::strcmp(*(atoms->atomname[i]), "CB") == 0) { atm.Cn[2] = i; } - else if ((strcmp(*(atoms->atomname[i]), "CG") == 0) || - (strcmp(*(atoms->atomname[i]), "CG1") == 0) || - (strcmp(*(atoms->atomname[i]), "OG") == 0) || - (strcmp(*(atoms->atomname[i]), "OG1") == 0) || - (strcmp(*(atoms->atomname[i]), "SG") == 0)) + else if ((std::strcmp(*(atoms->atomname[i]), "CG") == 0) || + (std::strcmp(*(atoms->atomname[i]), "CG1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OG") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OG1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "SG") == 0)) { atm.Cn[3] = i; } - else if ((strcmp(*(atoms->atomname[i]), "CD") == 0) || - (strcmp(*(atoms->atomname[i]), "CD1") == 0) || - (strcmp(*(atoms->atomname[i]), "SD") == 0) || - (strcmp(*(atoms->atomname[i]), "OD1") == 0) || - (strcmp(*(atoms->atomname[i]), "ND1") == 0)) + else if ((std::strcmp(*(atoms->atomname[i]), "CD") == 0) || + (std::strcmp(*(atoms->atomname[i]), "CD1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "SD") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OD1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "ND1") == 0)) { atm.Cn[4] = i; } /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */ - else if (bHChi && ((strcmp(*(atoms->atomname[i]), "HG") == 0) || - (strcmp(*(atoms->atomname[i]), "HG1") == 0)) ) + else if (bHChi && ((std::strcmp(*(atoms->atomname[i]), "HG") == 0) || + (std::strcmp(*(atoms->atomname[i]), "HG1") == 0)) ) { atm.Cn[4] = i; } - else if ((strcmp(*(atoms->atomname[i]), "CE") == 0) || - (strcmp(*(atoms->atomname[i]), "CE1") == 0) || - (strcmp(*(atoms->atomname[i]), "OE1") == 0) || - (strcmp(*(atoms->atomname[i]), "NE") == 0)) + else if ((std::strcmp(*(atoms->atomname[i]), "CE") == 0) || + (std::strcmp(*(atoms->atomname[i]), "CE1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "OE1") == 0) || + (std::strcmp(*(atoms->atomname[i]), "NE") == 0)) { atm.Cn[5] = i; } - else if ((strcmp(*(atoms->atomname[i]), "CZ") == 0) || - (strcmp(*(atoms->atomname[i]), "NZ") == 0)) + else if ((std::strcmp(*(atoms->atomname[i]), "CZ") == 0) || + (std::strcmp(*(atoms->atomname[i]), "NZ") == 0)) { atm.Cn[6] = i; } /* HChi flag here too */ - else if (bHChi && (strcmp(*(atoms->atomname[i]), "NH1") == 0)) + else if (bHChi && (std::strcmp(*(atoms->atomname[i]), "NH1") == 0)) { atm.Cn[7] = i; } @@ -153,13 +152,13 @@ t_dlist *mk_dlist(FILE *log, /* added by grs - special case for aromatics, whose chis above 2 are not real and produce rubbish output - so set back to -1 */ - if (strcmp(thisres, "PHE") == 0 || - strcmp(thisres, "TYR") == 0 || - strcmp(thisres, "PTR") == 0 || - strcmp(thisres, "TRP") == 0 || - strcmp(thisres, "HIS") == 0 || - strcmp(thisres, "HISA") == 0 || - strcmp(thisres, "HISB") == 0) + if (std::strcmp(thisres, "PHE") == 0 || + std::strcmp(thisres, "TYR") == 0 || + std::strcmp(thisres, "PTR") == 0 || + std::strcmp(thisres, "TRP") == 0 || + std::strcmp(thisres, "HIS") == 0 || + std::strcmp(thisres, "HISA") == 0 || + std::strcmp(thisres, "HISB") == 0) { for (ii = 5; ii <= 7; ii++) { @@ -423,7 +422,7 @@ int pr_trans(FILE *fp, int nl, t_dlist dl[], real dt, int Xi) fprintf(fp, "Residue\t&$\\chi_%d$\t\\\\\n", Xi+1); for (i = 0; (i < nl); i++) { - nn = dl[i].ntr[Xi]/dt; + nn = static_cast(dl[i].ntr[Xi]/dt); if (nn == 0) { diff --git a/src/gromacs/gmxana/edittop.c b/src/gromacs/gmxana/edittop.cpp similarity index 98% rename from src/gromacs/gmxana/edittop.c rename to src/gromacs/gmxana/edittop.cpp index 0338366db0..bc7db85fb7 100644 --- a/src/gromacs/gmxana/edittop.c +++ b/src/gromacs/gmxana/edittop.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -117,7 +117,7 @@ static void delete_from_interactions(t_idef *idef, int inr) static void delete_from_block(t_block *block, int inr) { /* Update block data structure */ - int i, i1, j1, j, k; + int i, i1, j; for (i = 0; (i < block->nr); i++) { @@ -193,8 +193,6 @@ static void delete_from_atoms(t_atoms *atoms, int inr) void delete_atom(t_topology *top, int inr) { - int k; - if ((inr < 0) || (inr >= top->atoms.nr)) { gmx_fatal(FARGS, "Delete_atom: inr (%d) not in %d .. %d", inr, 0, diff --git a/src/gromacs/gmxana/eigio.c b/src/gromacs/gmxana/eigio.cpp similarity index 100% rename from src/gromacs/gmxana/eigio.c rename to src/gromacs/gmxana/eigio.cpp diff --git a/src/gromacs/gmxana/eigio.h b/src/gromacs/gmxana/eigio.h index 7b2638e16d..82c3cb9c37 100644 --- a/src/gromacs/gmxana/eigio.h +++ b/src/gromacs/gmxana/eigio.h @@ -40,6 +40,11 @@ #include "gromacs/legacyheaders/typedefs.h" +#ifdef __cplusplus +extern "C" +{ +#endif + enum { eWXR_NO, eWXR_YES, eWXR_NOFIT }; @@ -81,4 +86,8 @@ int read_eigval (const char * fn, int eigvalnr[], real eigval[]); +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/gromacs/gmxana/fitahx.c b/src/gromacs/gmxana/fitahx.cpp similarity index 94% rename from src/gromacs/gmxana/fitahx.c rename to src/gromacs/gmxana/fitahx.cpp index 24b5a955d1..7616cc6523 100644 --- a/src/gromacs/gmxana/fitahx.c +++ b/src/gromacs/gmxana/fitahx.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -38,6 +38,8 @@ #include "fitahx.h" +#include + #include "gromacs/math/do_fit.h" #include "gromacs/math/vec.h" #include "gromacs/utility/fatalerror.h" @@ -97,8 +99,8 @@ real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[], for (i = 0; (i < nca); i++) { ai = caindex[i]; - xref[ai][XX] = rad*cos(phi0); - xref[ai][YY] = rad*sin(phi0); + xref[ai][XX] = rad*std::cos(phi0); + xref[ai][YY] = rad*std::sin(phi0); xref[ai][ZZ] = d*i; /* Set the mass to select that this Calpha contributes to fitting */ @@ -155,11 +157,11 @@ real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[], { rvec_sub(x[ai], xref[ai], dx); rms = iprod(dx, dx); - bb[i].rmsa += sqrt(rms); + bb[i].rmsa += std::sqrt(rms); bb[i].nrms++; trms += rms; mass[ai] = 0.0; } } - return sqrt(trms/nca); + return std::sqrt(trms/nca); } diff --git a/src/gromacs/gmxana/fitahx.h b/src/gromacs/gmxana/fitahx.h index 399028ef05..d7032a1640 100644 --- a/src/gromacs/gmxana/fitahx.h +++ b/src/gromacs/gmxana/fitahx.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -42,7 +42,16 @@ #include "gromacs/gmxana/hxprops.h" #include "gromacs/legacyheaders/typedefs.h" -extern real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[], - rvec x[], int nca, atom_id caindex[], gmx_bool bFit); +#ifdef __cplusplus +extern "C" +{ +#endif + +real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[], + rvec x[], int nca, atom_id caindex[], gmx_bool bFit); + +#ifdef __cplusplus +} +#endif #endif diff --git a/src/gromacs/gmxana/geminate.c b/src/gromacs/gmxana/geminate.c deleted file mode 100644 index d06a4b41d8..0000000000 --- a/src/gromacs/gmxana/geminate.c +++ /dev/null @@ -1,767 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public License - * as published by the Free Software Foundation; either version 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "geminate.h" - -#include -#include - -#include "gromacs/legacyheaders/typedefs.h" -#include "gromacs/math/vec.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/gmxomp.h" -#include "gromacs/utility/smalloc.h" - -static void missing_code_message() -{ - fprintf(stderr, "You have requested code to run that is deprecated.\n"); - fprintf(stderr, "Revert to an older GROMACS version or help in porting the code.\n"); -} - -/* The first few sections of this file contain functions that were adopted, - * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se, - * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url). - * This is also the case with the function eq10v2(). - * - * The parts menetioned in the previous paragraph were contributed under the BSD license. - */ - - -/* This first part is from complex.c which I recieved from Omer Markowitch. - * - Erik Marklund - * - * ------------- from complex.c ------------- */ - -/* Complexation of a paired number (r,i) */ -static gem_complex gem_cmplx(double r, double i) -{ - gem_complex value; - value.r = r; - value.i = i; - return value; -} - -/* Complexation of a real number, x */ -static gem_complex gem_c(double x) -{ - gem_complex value; - value.r = x; - value.i = 0; - return value; -} - -/* Magnitude of a complex number z */ -static double gem_cx_abs(gem_complex z) -{ - return (sqrt(z.r*z.r+z.i*z.i)); -} - -/* Addition of two complex numbers z1 and z2 */ -static gem_complex gem_cxadd(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r+z2.r; - value.i = z1.i+z2.i; - return value; -} - -/* Addition of a complex number z1 and a real number r */ -static gem_complex gem_cxradd(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r + r; - value.i = z.i; - return value; -} - -/* Subtraction of two complex numbers z1 and z2 */ -static gem_complex gem_cxsub(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r-z2.r; - value.i = z1.i-z2.i; - return value; -} - -/* Multiplication of two complex numbers z1 and z2 */ -static gem_complex gem_cxmul(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r*z2.r-z1.i*z2.i; - value.i = z1.r*z2.i+z1.i*z2.r; - return value; -} - -/* Square of a complex number z */ -static gem_complex gem_cxsq(gem_complex z) -{ - gem_complex value; - value.r = z.r*z.r-z.i*z.i; - value.i = z.r*z.i*2.; - return value; -} - -/* multiplication of a complex number z and a real number r */ -static gem_complex gem_cxrmul(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r*r; - value.i = z.i*r; - return value; -} - -/* Division of two complex numbers z1 and z2 */ -static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2) -{ - gem_complex value; - double num; - num = z2.r*z2.r+z2.i*z2.i; - if (num == 0.) - { - fprintf(stderr, "ERROR in gem_cxdiv function\n"); - } - value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num; - return value; -} - -/* division of a complex z number by a real number x */ -static gem_complex gem_cxrdiv(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r/r; - value.i = z.i/r; - return value; -} - -/* division of a real number r by a complex number x */ -static gem_complex gem_rxcdiv(double r, gem_complex z) -{ - gem_complex value; - double f; - f = r/(z.r*z.r+z.i*z.i); - value.r = f*z.r; - value.i = -f*z.i; - return value; -} - -/* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/ -static gem_complex gem_cxdexp(gem_complex z) -{ - gem_complex value; - double exp_z_r; - exp_z_r = exp(z.r); - value.r = exp_z_r*cos(z.i); - value.i = exp_z_r*sin(z.i); - return value; -} - -/* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z), */ -/* where -PI < Arg(z) < PI */ -static gem_complex gem_cxlog(gem_complex z) -{ - gem_complex value; - double mag2; - mag2 = z.r*z.r+z.i*z.i; - if (mag2 < 0.) - { - fprintf(stderr, "ERROR in gem_cxlog func\n"); - } - value.r = log(sqrt(mag2)); - if (z.r == 0.) - { - value.i = PI/2.; - if (z.i < 0.) - { - value.i = -value.i; - } - } - else - { - value.i = atan2(z.i, z.r); - } - return value; -} - -/* Square root of a complex number z = |z| exp(I*the) -- z^(1/2) */ -/* z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */ -/* where 0 < the < 2*PI */ -static gem_complex gem_cxdsqrt(gem_complex z) -{ - gem_complex value; - double sq; - sq = gem_cx_abs(z); - value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */ - value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */ - if (z.i < 0.) - { - value.r = -value.r; - } - return value; -} - -/* Complex power of a complex number z1^z2 */ -static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2)); - return value; -} - -/* ------------ end of complex.c ------------ */ - -/* This next part was derived from cubic.c, also received from Omer Markovitch. - * ------------- from cubic.c ------------- */ - -/* Solver for a cubic equation: x^3-a*x^2+b*x-c=0 */ -static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam, - double a, double b, double c) -{ - double t1, t2, two_3, temp; - gem_complex ctemp, ct3; - - two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c; - temp = 4.*t1*t1*t1+t2*t2; - - ctemp = gem_cmplx(temp, 0.); ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp)); - ct3 = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.)); - - ctemp = gem_rxcdiv(-two_3*t1/3., ct3); - ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3)); - - *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp); - - ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.)))); - ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam)); - ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam); - *al = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5); - *be = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5); -} - -/* ------------ end of cubic.c ------------ */ - -/* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch. - * ------------- from [cr]error.c ------------- */ - -/************************************************************/ -/* Real valued error function and related functions */ -/* x, y : real variables */ -/* erf(x) : error function */ -/* erfc(x) : complementary error function */ -/* omega(x) : exp(x*x)*erfc(x) */ -/* W(x,y) : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */ -/************************************************************/ - -/*---------------------------------------------------------------------------*/ -/* Utilzed the series approximation of erf(x) */ -/* Relative error=|err(x)|/erf(x)=6 series sum deteriorates -> Asymptotic series used instead */ -/*---------------------------------------------------------------------------*/ -/* This gives erfc(z) correctly only upto >10-15 */ - -static double gem_erf(double x) -{ - double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12; - temp = x; - sum = temp; - xm = 26.; - x2 = x*x; - x4 = x2*x2; - x6 = x4*x2; - x8 = x6*x2; - x10 = x8*x2; - x12 = x10*x2; - exp2 = exp(-x2); - if (x <= xm) - { - for (n = 1.; n <= 2000.; n += 1.) - { - temp *= 2.*x2/(2.*n+1.); - sum += temp; - if (fabs(temp/sum) < 1.E-16) - { - break; - } - } - - if (n >= 2000.) - { - fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n); - } - sum *= 2./sPI*exp2; - } - else - { - /* from the asymptotic expansion of experfc(x) */ - sum = (1. - 0.5/x2 + 0.75/x4 - - 1.875/x6 + 6.5625/x8 - - 29.53125/x10 + 162.421875/x12) - / sPI/x; - sum *= exp2; /* now sum is erfc(x) */ - sum = -sum+1.; - } - return x >= 0.0 ? sum : -sum; -} - -/* Result --> Alex's code for erfc and experfc behaves better than this */ -/* complementray error function. Returns 1.-erf(x) */ -static double gem_erfc(double x) -{ - double t, z, ans; - z = fabs(x); - t = 1.0/(1.0+0.5*z); - - ans = t * exp(-z*z - 1.26551223 + - t*(1.00002368 + - t*(0.37409196 + - t*(0.09678418 + - t*(-0.18628806 + - t*(0.27886807 + - t*(-1.13520398 + - t*(1.48851587 + - t*(-0.82215223 + - t*0.17087277))))))))); - - return x >= 0.0 ? ans : 2.0-ans; -} - -/* omega(x)=exp(x^2)erfc(x) */ -static double gem_omega(double x) -{ - double xm, ans, xx, x4, x6, x8, x10, x12; - xm = 26; - xx = x*x; - x4 = xx*xx; - x6 = x4*xx; - x8 = x6*xx; - x10 = x8*xx; - x12 = x10*xx; - - if (x <= xm) - { - ans = exp(xx)*gem_erfc(x); - } - else - { - /* Asymptotic expansion */ - ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x; - } - return ans; -} - -/*---------------------------------------------------------------------------*/ -/* Utilzed the series approximation of erf(z=x+iy) */ -/* Relative error=|err(z)|/|erf(z)|kD; - D = params->D; - r = params->sigma; - a = (1.0 + ka/kD) * sqrt(D)/r; - b = kd; - c = kd * sqrt(D)/r; - - gem_solve(&alpha, &beta, &gamma, a, b, c); - /* Finding the 3 roots */ - - sumimaginary = 0; - part1 = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta, gamma), gem_cxsub(beta, gamma))); /* 1(2+3)(2-3) */ - part2 = gem_cxmul(beta, gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha))); /* 2(3+1)(3-1) */ - part3 = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta))); /* 3(1+2)(1-2) */ - part4 = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */ - -#pragma omp parallel for \ - private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \ - reduction(+:sumimaginary) \ - default(shared) \ - schedule(guided) - for (i = 0; i < manytimes; i++) - { - tsqrt = sqrt(time[i]); - oma = gem_comega(gem_cxrmul(alpha, tsqrt)); - c1 = gem_cxmul(oma, gem_cxdiv(part1, part4)); - omb = gem_comega(gem_cxrmul(beta, tsqrt)); - c2 = gem_cxmul(omb, gem_cxdiv(part2, part4)); - omc = gem_comega(gem_cxrmul(gamma, tsqrt)); - c3 = gem_cxmul(omc, gem_cxdiv(part3, part4)); - c4.r = c1.r+c2.r+c3.r; - c4.i = c1.i+c2.i+c3.i; - theoryCt[i] = c4.r; - sumimaginary += c4.i * c4.i; - } - - return sumimaginary; - -} /* eq10v2 */ - -/* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */ -static double getLogIndex(const int i, const t_gemParams *params) -{ - return gmx_expm1(((double)(i)) * params->logQuota); -} - -extern t_gemParams *init_gemParams(const double sigma, const double D, - const real *t, const int len, const int nFitPoints, - const real begFit, const real endFit, - const real ballistic, const int nBalExp) -{ - double tDelta; - t_gemParams *p; - - snew(p, 1); - - /* A few hardcoded things here. For now anyway. */ -/* p->ka_min = 0; */ -/* p->ka_max = 100; */ -/* p->dka = 10; */ -/* p->kd_min = 0; */ -/* p->kd_max = 2; */ -/* p->dkd = 0.2; */ - p->ka = 0; - p->kd = 0; -/* p->lsq = -1; */ -/* p->lifetime = 0; */ - p->sigma = sigma*10; /* Omer uses Å, not nm */ -/* p->lsq_old = 99999; */ - p->D = sqcm_per_s_to_sqA_per_ps(D); - p->kD = 4*acos(-1.0)*sigma*p->D; - - - /* Parameters used by calcsquare(). Better to calculate them - * here than in calcsquare every time it's called. */ - p->len = len; -/* p->logAfterTime = logAfterTime; */ - tDelta = (t[len-1]-t[0]) / len; - if (tDelta <= 0) - { - gmx_fatal(FARGS, "Time between frames is non-positive!"); - } - else - { - p->tDelta = tDelta; - } - - p->nExpFit = nBalExp; -/* p->nLin = logAfterTime / tDelta; */ - p->nFitPoints = nFitPoints; - p->begFit = begFit; - p->endFit = endFit; - p->logQuota = (double)(log(p->len))/(p->nFitPoints-1); -/* if (p->nLin <= 0) { */ -/* fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */ -/* sfree(p); */ -/* return NULL; */ -/* } */ -/* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */ -/* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */ -/* p->logPF = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */ -/* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */ - -/* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */ - p->ballistic = ballistic; - return p; -} - -/* There was a misunderstanding regarding the fitting. From our - * recent correspondence it appears that Omer's code require - * the ACF data on a log-scale and does not operate on the raw data. - * This needs to be redone in gemFunc_residual() as well as in the - * t_gemParams structure. */ - -static real* d2r(const double *d, const int nn); - -extern real fitGemRecomb(double gmx_unused *ct, - double gmx_unused *time, - double gmx_unused **ctFit, - const int gmx_unused nData, - t_gemParams gmx_unused *params) -{ - - int nThreads, i, iter, status, maxiter; - real size, d2, tol, *dumpdata; - size_t p, n; - gemFitData *GD; - char *dumpstr, dumpname[128]; - - missing_code_message(); - return -1; - -} - - -/* Removes the ballistic term from the beginning of the ACF, - * just like in Omer's paper. - */ -extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative) -{ - - /* Fit with 4 exponentials and one constant term, - * subtract the fatest exponential. */ - - int nData, i, status, iter; - balData *BD; - double *guess, /* Initial guess. */ - *A, /* The fitted parameters. (A1, B1, A2, B2,... C) */ - a[2], - ddt[2]; - gmx_bool sorted; - size_t n; - size_t p; - - nData = 0; - do - { - nData++; - } - while (t[nData] < tMax+t[0] && nData < len); - - p = nexp*2+1; /* Number of parameters. */ - - missing_code_message(); - return; -} - -extern void dumpN(const real *e, const int nn, char *fn) -{ - /* For debugging only */ - int i; - FILE *f; - char standardName[] = "Nt.xvg"; - if (fn == NULL) - { - fn = standardName; - } - - f = fopen(fn, "w"); - fprintf(f, - "@ type XY\n" - "@ xaxis label \"Frame\"\n" - "@ yaxis label \"N\"\n" - "@ s0 line type 3\n"); - - for (i = 0; i < nn; i++) - { - fprintf(f, "%-10i %-g\n", i, e[i]); - } - - fclose(f); -} - -static real* d2r(const double *d, const int nn) -{ - real *r; - int i; - - snew(r, nn); - for (i = 0; i < nn; i++) - { - r[i] = (real)d[i]; - } - - return r; -} - -static void _patchBad(double *ct, int n, double dy) -{ - /* Just do lin. interpolation for now. */ - int i; - - for (i = 1; i < n; i++) - { - ct[i] = ct[0]+i*dy; - } -} - -static void patchBadPart(double *ct, int n) -{ - _patchBad(ct, n, (ct[n] - ct[0])/n); -} - -static void patchBadTail(double *ct, int n) -{ - _patchBad(ct+1, n-1, ct[1]-ct[0]); - -} - -extern void fixGemACF(double *ct, int len) -{ - int i, j, b, e; - gmx_bool bBad; - - /* Let's separate two things: - * - identification of bad parts - * - patching of bad parts. - */ - - b = 0; /* Start of a bad stretch */ - e = 0; /* End of a bad stretch */ - bBad = FALSE; - - /* An acf of binary data must be one at t=0. */ - if (fabs(ct[0]-1.0) > 1e-6) - { - ct[0] = 1.0; - fprintf(stderr, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct[0]-1.0)); - } - - for (i = 0; i < len; i++) - { - -#ifdef HAS_ISFINITE - if (isfinite(ct[i])) -#elif defined(HAS__ISFINITE) - if (_isfinite(ct[i])) -#else - if (1) -#endif - { - if (!bBad) - { - /* Still on a good stretch. Proceed.*/ - continue; - } - - /* Patch up preceding bad stretch. */ - if (i == (len-1)) - { - /* It's the tail */ - if (b <= 1) - { - gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting."); - } - patchBadTail(&(ct[b-2]), (len-b)+1); - } - - e = i; - patchBadPart(&(ct[b-1]), (e-b)+1); - bBad = FALSE; - } - else - { - if (!bBad) - { - b = i; - - bBad = TRUE; - } - } - } -} diff --git a/src/gromacs/gmxana/geminate.h b/src/gromacs/gmxana/geminate.h deleted file mode 100644 index 91bb097421..0000000000 --- a/src/gromacs/gmxana/geminate.h +++ /dev/null @@ -1,184 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public License - * as published by the Free Software Foundation; either version 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#ifndef _GEMINATE_H -#define _GEMINATE_H - -#include - -#include "gromacs/utility/basedefinitions.h" -#include "gromacs/utility/real.h" - -enum { - gemNULL, gemNONE, gemDD, gemAD, gemAA, gemA4, gemNR -}; -static const char *gemType[] = {NULL, "none", "dd", "ad", "aa", "a4", NULL}; - -/* The first few sections of this file contain functions that were adopted, - * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se, - * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url). - * This is also the case with the function eq10v2() in geminate.c. - * - * The parts menetioned in the previous paragraph were contributed under a BSD license. - */ - -/* This first part is derived from complex.c which I recieved from Omer Markowitch. - * - Erik Marklund - * - * ------------- from complex.c ------------- */ - -#include -/* definition of PI */ -#ifndef PI -#define PI (acos(-1.0)) -#endif - -/* definition of the type complex */ -typedef struct -{ - double r, i; -} gem_complex; - - -/* ------------ end of complex.c ------------ */ - -/* This next part was derived from cerror.c and rerror.c, - * also received from Omer Markovitch. - * ------------- from [cr]error.c ------------- */ - -#ifndef sPI -#define sPI (sqrt(PI)) -#endif - -/* ------------ end of [cr]error.c ------------ */ - -/* ///////////////// REVERSIBLE GEMINATE RECOMBINATION /////////////////// - * Here follow routines and structs for reversible geminate recombination. - */ - -typedef struct { - size_t n; - double *y; - double tDelta; - int nexp; -} balData; - - -typedef struct { - /* Used in the rewritten version of Omer's gem. recomb. analysis */ - double ka, kd; /* -- || -- results[] */ - double sigma; /* -- || -- sigma */ - double D; /* The diffusion coefficient */ - double kD; /* Replaces kD in analyse_corr_gem3d() */ - - /* The following members are for calcsquare() and takeAwayBallistic() */ - double tDelta; /* Time between frames */ - /* double logAfterTime; /\* Time after which we do the lsq calculations on a logarithmic timescale. *\/ */ - int nFitPoints; /* Number of points to take from the ACF for fitting */ - double begFit; /* Fit from this time (ps) */ - double endFit; /* Fit up to this time (ps) */ -/* double logDelta; */ -/* double logPF; */ -/* To get an equal number of points in the lin and log regime, - * we'll use logDelta to calculate where to sample the ACF. - * if i and j are indices in the linear and log regime, then: - * j = Cexp(A(i+nLin)), - * where C = (nLin**2 / len) and A = log(len/nLin) / nLin. - * This expands to j = (nLin**2 / len) * exp((i+nLin) * log(len/nLin) / nLin). - * We'll save part of our brains and some computing time if we pre-calculate - * 1) logDelta = log(len/nLin) / nLin - * 2) logPF = nLin**2 / len - * and get j = logPF * exp((i+nLin) * logDelta). */ - - /* I will redo this for a fit done entirely in log-log. - * j' = j+1 - * nFitPoints' = nFitPoints-1 - * - * j' = Cexp(Ai) - * (j'= 1 <=> i=0) - * => C=1 - * (j'=len <=> i=nFitPoints') - * => A=log(len)/nFitPoints' - * => j = exp(i*log(len)/(nFitPoints-1)) -1 - **/ -/* #define GETLOGINDEX(i,params) (params)->logPF * exp(((i)+(params)->nLin) * (params)->logDelta) - */ - double logQuota; - int nLin; /* Number of timepoints in the linear regime */ - int len; /* Length of time and ct arrays */ - int nExpFit; /* Number of exponentials to fit */ - real ballistic; /* Time before which the ballistic term should be fitted */ - gmx_bool bDt; /* TRUE => use time derivative at time 0 - * to find fastest component. - * FALSE => use coefficient in exponenetial - * to find fastest component. */ -} t_gemParams; - - -typedef struct { - size_t n; /* Number of data points (lin-log) */ - double *y; /* The datapoints */ - double *ctTheory; /* The theoretical ct which will be built by gemFunc_f. */ - double *LinLog; - double *time; - double ka; - double kd; - double tDelta; /* time difference between subsequent datapoints */ - size_t nData; /* real size of the data */ - int *logtime; - double *doubleLogTime; - t_gemParams *params; -} gemFitData; - -extern void takeAwayBallistic(double *ct, double *t, - int len, real tMax, - int nexp, gmx_bool bDerivative); - - -extern t_gemParams *init_gemParams(const double sigma, const double D, - const real *t, const int len, const int nFitPoints, - const real begFit, const real endFit, - const real ballistic, const int nBalExp); - -/* Fit to geminate recombination model. - Returns root mean square error of fit. */ -extern real fitGemRecomb(double *ct, double *time, double **ctFit, - const int nData, t_gemParams *params); - -extern void dumpN(const real *e, const int nn, char *fn); - -/* Fix NaN that might appear in the theoretical acf. */ -extern void fixGemACF(double *ct, int len); - -#endif diff --git a/src/gromacs/gmxana/gmx_analyze.c b/src/gromacs/gmxana/gmx_analyze.c index 93ccf37dd6..5d1c8c4a59 100644 --- a/src/gromacs/gmxana/gmx_analyze.c +++ b/src/gromacs/gmxana/gmx_analyze.c @@ -45,7 +45,6 @@ #include "gromacs/correlationfunctions/expfit.h" #include "gromacs/correlationfunctions/integrate.h" #include "gromacs/fileio/xvgr.h" -#include "gromacs/gmxana/geminate.h" #include "gromacs/gmxana/gmx_ana.h" #include "gromacs/gmxana/gstat.h" #include "gromacs/legacyheaders/copyrite.h" @@ -904,122 +903,6 @@ static void do_fit(FILE *out, int n, gmx_bool bYdy, } } -static void do_ballistic(const char *balFile, int nData, - real *t, real **val, int nSet, - real balTime, int nBalExp, - const output_env_t oenv) -{ - double **ctd = NULL, *td = NULL; - t_gemParams *GP = init_gemParams(0, 0, t, nData, 0, 0, 0, balTime, nBalExp); - static char *leg[] = {"Ac'(t)"}; - FILE *fp; - int i, set; - - if (GP->ballistic/GP->tDelta >= GP->nExpFit*2+1) - { - snew(ctd, nSet); - snew(td, nData); - - fp = xvgropen(balFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv); - xvgr_legend(fp, asize(leg), (const char**)leg, oenv); - - for (set = 0; set < nSet; set++) - { - snew(ctd[set], nData); - for (i = 0; i < nData; i++) - { - ctd[set][i] = (double)val[set][i]; - if (set == 0) - { - td[i] = (double)t[i]; - } - } - - takeAwayBallistic(ctd[set], td, nData, GP->ballistic, GP->nExpFit, GP->bDt); - } - - for (i = 0; i < nData; i++) - { - fprintf(fp, " %g", t[i]); - for (set = 0; set < nSet; set++) - { - fprintf(fp, " %g", ctd[set][i]); - } - fprintf(fp, "\n"); - } - - - for (set = 0; set < nSet; set++) - { - sfree(ctd[set]); - } - sfree(ctd); - sfree(td); - xvgrclose(fp); - } - else - { - printf("Number of data points is less than the number of parameters to fit\n." - "The system is underdetermined, hence no ballistic term can be found.\n\n"); - } -} - -static void do_geminate(const char *gemFile, int nData, - real *t, real **val, int nSet, - const real D, const real rcut, const real balTime, - const int nFitPoints, const real begFit, const real endFit, - const output_env_t oenv) -{ - double **ctd = NULL, **ctdGem = NULL, *td = NULL; - t_gemParams *GP = init_gemParams(rcut, D, t, nData, nFitPoints, - begFit, endFit, balTime, 1); - const char *leg[] = {"Ac\\sgem\\N(t)"}; - FILE *fp; - int i, set; - - snew(ctd, nSet); - snew(ctdGem, nSet); - snew(td, nData); - - fp = xvgropen(gemFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv); - xvgr_legend(fp, asize(leg), leg, oenv); - - for (set = 0; set < nSet; set++) - { - snew(ctd[set], nData); - snew(ctdGem[set], nData); - for (i = 0; i < nData; i++) - { - ctd[set][i] = (double)val[set][i]; - if (set == 0) - { - td[i] = (double)t[i]; - } - } - fitGemRecomb(ctd[set], td, &(ctd[set]), nData, GP); - } - - for (i = 0; i < nData; i++) - { - fprintf(fp, " %g", t[i]); - for (set = 0; set < nSet; set++) - { - fprintf(fp, " %g", ctdGem[set][i]); - } - fprintf(fp, "\n"); - } - - for (set = 0; set < nSet; set++) - { - sfree(ctd[set]); - sfree(ctdGem[set]); - } - sfree(ctd); - sfree(ctdGem); - sfree(td); - xvgrclose(fp); -} - static void print_fitted_function(const char *fitfile, const char *fn_fitted, gmx_bool bXYdy, @@ -1129,20 +1012,6 @@ int gmx_analyze(int argc, char *argv[]) "The complete derivation is given in", "B. Hess, J. Chem. Phys. 116:209-217, 2002.[PAR]", - "Option [TT]-bal[tt] finds and subtracts the ultrafast \"ballistic\"", - "component from a hydrogen bond autocorrelation function by the fitting", - "of a sum of exponentials, as described in e.g.", - "O. Markovitch, J. Chem. Phys. 129:084505, 2008. The fastest term", - "is the one with the most negative coefficient in the exponential,", - "or with [TT]-d[tt], the one with most negative time derivative at time 0.", - "[TT]-nbalexp[tt] sets the number of exponentials to fit.[PAR]", - - "Option [TT]-gem[tt] fits bimolecular rate constants ka and kb", - "(and optionally kD) to the hydrogen bond autocorrelation function", - "according to the reversible geminate recombination model. Removal of", - "the ballistic component first is strongly advised. The model is presented in", - "O. Markovitch, J. Chem. Phys. 129:084505, 2008.[PAR]", - "Option [TT]-filter[tt] prints the RMS high-frequency fluctuation", "of each set and over all sets with respect to a filtered average.", "The filter is proportional to cos([GRK]pi[grk] t/len) where t goes from -len/2", @@ -1171,8 +1040,8 @@ int gmx_analyze(int argc, char *argv[]) static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE; static gmx_bool bEESEF = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE; static gmx_bool bIntegrate = FALSE, bRegression = FALSE, bLuzar = FALSE, bLuzarError = FALSE; - static int nsets_in = 1, d = 1, nb_min = 4, resol = 10, nBalExp = 4, nFitPoints = 100; - static real temp = 298.15, fit_start = 1, fit_end = 60, balTime = 0.2, diffusion = 5e-5, rcut = 0.35; + static int nsets_in = 1, d = 1, nb_min = 4, resol = 10, nFitPoints = 100; + static real temp = 298.15, fit_start = 1, fit_end = 60, diffusion = 5e-5, rcut = 0.35; /* must correspond to enum avbar* declared at beginning of file */ static const char *avbar_opt[avbarNR+1] = { @@ -1234,18 +1103,6 @@ int gmx_analyze(int argc, char *argv[]) "Subtract the average before autocorrelating" }, { "-oneacf", FALSE, etBOOL, {&bAverCorr}, "Calculate one ACF over all sets" }, - { "-nbalexp", FALSE, etINT, {&nBalExp}, - "HIDDENNumber of exponentials to fit to the ultrafast component" }, - { "-baltime", FALSE, etREAL, {&balTime}, - "HIDDENTime up to which the ballistic component will be fitted" }, -/* { "-gemnp", FALSE, etINT, {&nFitPoints}, */ -/* "HIDDENNumber of data points taken from the ACF to use for fitting to rev. gem. recomb. model."}, */ -/* { "-rcut", FALSE, etREAL, {&rcut}, */ -/* "Cut-off for hydrogen bonds in geminate algorithms" }, */ -/* { "-gemtype", FALSE, etENUM, {gemType}, */ -/* "What type of gminate recombination to use"}, */ -/* { "-D", FALSE, etREAL, {&diffusion}, */ -/* "The self diffusion coefficient which is used for the reversible geminate recombination model."} */ }; #define NPA asize(pa) @@ -1253,7 +1110,7 @@ int gmx_analyze(int argc, char *argv[]) int n, nlast, s, nset, i, j = 0; real **val, *t, dt, tot, error; double *av, *sig, cum1, cum2, cum3, cum4, db; - const char *acfile, *msdfile, *ccfile, *distfile, *avfile, *eefile, *balfile, *gemfile, *fitfile; + const char *acfile, *msdfile, *ccfile, *distfile, *avfile, *eefile, *fitfile; output_env_t oenv; t_filenm fnm[] = { @@ -1264,9 +1121,7 @@ int gmx_analyze(int argc, char *argv[]) { efXVG, "-dist", "distr", ffOPTWR }, { efXVG, "-av", "average", ffOPTWR }, { efXVG, "-ee", "errest", ffOPTWR }, - { efXVG, "-bal", "ballisitc", ffOPTWR }, { efXVG, "-fitted", "fitted", ffOPTWR }, -/* { efXVG, "-gem", "geminate", ffOPTWR }, */ { efLOG, "-g", "fitlog", ffOPTWR } }; #define NFILE asize(fnm) @@ -1289,11 +1144,6 @@ int gmx_analyze(int argc, char *argv[]) distfile = opt2fn_null("-dist", NFILE, fnm); avfile = opt2fn_null("-av", NFILE, fnm); eefile = opt2fn_null("-ee", NFILE, fnm); - balfile = opt2fn_null("-bal", NFILE, fnm); -/* gemfile = opt2fn_null("-gem",NFILE,fnm); */ - /* When doing autocorrelation we don't want a fitlog for fitting - * the function itself (not the acf) when the user did not ask for it. - */ if (opt2parg_bSet("-fitfn", npargs, ppa) && acfile == NULL) { fitfile = opt2fn("-g", NFILE, fnm); @@ -1450,13 +1300,6 @@ int gmx_analyze(int argc, char *argv[]) estimate_error(eefile, nb_min, resol, n, nset, av, sig, val, dt, bEeFitAc, bEESEF, bEENLC, oenv); } - if (balfile) - { - do_ballistic(balfile, n, t, val, nset, balTime, nBalExp, oenv); - } -/* if (gemfile) */ -/* do_geminate(gemfile,n,t,val,nset,diffusion,rcut,balTime, */ -/* nFitPoints, fit_start, fit_end, oenv); */ if (bPower) { power_fit(n, nset, val, t); diff --git a/src/gromacs/gmxana/gmx_hbond.c b/src/gromacs/gmxana/gmx_hbond.c index 73b17d1a05..8bde9f1ac0 100644 --- a/src/gromacs/gmxana/gmx_hbond.c +++ b/src/gromacs/gmxana/gmx_hbond.c @@ -49,7 +49,6 @@ #include "gromacs/fileio/tpxio.h" #include "gromacs/fileio/trxio.h" #include "gromacs/fileio/xvgr.h" -#include "gromacs/gmxana/geminate.h" #include "gromacs/gmxana/gmx_ana.h" #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/legacyheaders/macros.h" @@ -66,12 +65,6 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/snprintf.h" -#include "geminate.h" - -/*#define HAVE_NN_LOOPS*/ - -typedef short int t_E; -typedef int t_EEst; #define max_hx 7 typedef int t_hx[max_hx]; #define NRHXTYPES max_hx @@ -96,9 +89,6 @@ enum { enum { noDA, ACC, DON, DA, INGROUP }; -enum { - NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR -}; static const char *grpnames[grNR] = {"0", "1", "I" }; @@ -168,39 +158,8 @@ typedef struct { h_id *nhbonds; /* The number of HBs per H at current */ } t_donors; -/* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/ -#define PSTYPE int - -typedef struct { - int len; /* The length of frame and p. */ - int *frame; /* The frames at which transitio*/ - PSTYPE *p; -} t_pShift; - -typedef struct { - /* Periodicity history. Used for the reversible geminate recombination. */ - t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap: - * pHist[d][a]. We can safely assume that the same - * periodic shift holds for all hydrogens of a da-pair. - * - * Nowadays it only stores TRANSITIONS, and not the shift at every frame. - * That saves a LOT of memory, an hopefully kills a mysterious bug where - * pHist gets contaminated. */ - - PSTYPE nper; /* The length of p2i */ - ivec *p2i; /* Maps integer to periodic shift for a pair.*/ - matrix P; /* Projection matrix to find the box shifts. */ - int gemtype; /* enumerated type */ -} t_gemPeriod; - -typedef struct { - int nframes; - int *Etot; /* Total energy for each frame */ - t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */ -} t_hbEmap; - typedef struct { - gmx_bool bHBmap, bDAnr, bGem; + gmx_bool bHBmap, bDAnr; int wordlen; /* The following arrays are nframes long */ int nframes, max_frames, maxhydro; @@ -215,128 +174,15 @@ typedef struct { /* This holds a matrix with all possible hydrogen bonds */ int nrhb, nrdist; t_hbond ***hbmap; -#ifdef HAVE_NN_LOOPS - t_hbEmap hbE; -#endif - /* For parallelization reasons this will have to be a pointer. - * Otherwise discrepancies may arise between the periodicity data - * seen by different threads. */ - t_gemPeriod *per; } t_hbdata; -static void clearPshift(t_pShift *pShift) -{ - if (pShift->len > 0) - { - sfree(pShift->p); - sfree(pShift->frame); - pShift->len = 0; - } -} - -static void calcBoxProjection(matrix B, matrix P) -{ - const int vp[] = {XX, YY, ZZ}; - int i, j; - int m, n; - matrix M, N, U; - - for (i = 0; i < 3; i++) - { - m = vp[i]; - for (j = 0; j < 3; j++) - { - n = vp[j]; - U[m][n] = i == j ? 1 : 0; - } - } - m_inv(B, M); - for (i = 0; i < 3; i++) - { - m = vp[i]; - mvmul(M, U[m], P[m]); - } - transpose(P, N); -} - -static void calcBoxDistance(matrix P, rvec d, ivec ibd) -{ - /* returns integer distance in box coordinates. - * P is the projection matrix from cartesian coordinates - * obtained with calcBoxProjection(). */ - int i; - rvec bd; - mvmul(P, d, bd); - /* extend it by 0.5 in all directions since (int) rounds toward 0.*/ - for (i = 0; i < 3; i++) - { - bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5); - } - ibd[XX] = (int)bd[XX]; - ibd[YY] = (int)bd[YY]; - ibd[ZZ] = (int)bd[ZZ]; -} - /* Changed argument 'bMerge' into 'oneHB' below, * since -contact should cause maxhydro to be 1, * not just -merge. * - Erik Marklund May 29, 2006 */ -static PSTYPE periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) -{ - /* Try to merge hbonds on the fly. That means that if the - * acceptor and donor are mergable, then: - * 1) store the hb-info so that acceptor id > donor id, - * 2) add the periodic shift in pairs, so that [-x,-y,-z] is - * stored in per.p2i[] whenever acceptor id < donor id. - * Note that [0,0,0] should already be the first element of per.p2i - * by the time this function is called. */ - - /* daSwap is TRUE if the donor and acceptor were swapped. - * If so, then the negative vector should be used. */ - PSTYPE i; - - if (per->p2i == NULL || per->nper == 0) - { - gmx_fatal(FARGS, "'per' not initialized properly."); - } - for (i = 0; i < per->nper; i++) - { - if (r[XX] == per->p2i[i][XX] && - r[YY] == per->p2i[i][YY] && - r[ZZ] == per->p2i[i][ZZ]) - { - return i; - } - } - /* Not found apparently. Add it to the list! */ - /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */ - -#pragma omp critical - { - if (!per->p2i) - { - fprintf(stderr, "p2i not initialized. This shouldn't happen!\n"); - snew(per->p2i, 1); - } - else - { - srenew(per->p2i, per->nper+2); - } - copy_ivec(r, per->p2i[per->nper]); - (per->nper)++; - - /* Add the mirror too. It's rather likely that it'll be needed. */ - per->p2i[per->nper][XX] = -r[XX]; - per->p2i[per->nper][YY] = -r[YY]; - per->p2i[per->nper][ZZ] = -r[ZZ]; - (per->nper)++; - } /* omp critical */ - return per->nper - 1 - (daSwap ? 0 : 1); -} - -static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode) +static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB) { t_hbdata *hb; @@ -344,7 +190,6 @@ static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_ hb->wordlen = 8*sizeof(unsigned int); hb->bHBmap = bHBmap; hb->bDAnr = bDAnr; - hb->bGem = bGem; if (oneHB) { hb->maxhydro = 1; @@ -353,9 +198,6 @@ static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_ { hb->maxhydro = MAXHYDRO; } - snew(hb->per, 1); - hb->per->gemtype = bGem ? gemmode : 0; - return hb; } @@ -378,295 +220,6 @@ static void mk_hbmap(t_hbdata *hb) } } -/* Consider redoing pHist so that is only stores transitions between - * periodicities and not the periodicity for all frames. This eats heaps of memory. */ -static void mk_per(t_hbdata *hb) -{ - int i, j; - if (hb->bGem) - { - snew(hb->per->pHist, hb->d.nrd); - for (i = 0; i < hb->d.nrd; i++) - { - snew(hb->per->pHist[i], hb->a.nra); - if (hb->per->pHist[i] == NULL) - { - gmx_fatal(FARGS, "Could not allocate enough memory for per->pHist"); - } - for (j = 0; j < hb->a.nra; j++) - { - clearPshift(&(hb->per->pHist[i][j])); - } - } - /* add the [0,0,0] shift to element 0 of p2i. */ - snew(hb->per->p2i, 1); - clear_ivec(hb->per->p2i[0]); - hb->per->nper = 1; - } -} - -#ifdef HAVE_NN_LOOPS -static void mk_hbEmap (t_hbdata *hb, int n0) -{ - int i, j, k; - hb->hbE.E = NULL; - hb->hbE.nframes = 0; - snew(hb->hbE.E, hb->d.nrd); - for (i = 0; i < hb->d.nrd; i++) - { - snew(hb->hbE.E[i], hb->a.nra); - for (j = 0; j < hb->a.nra; j++) - { - snew(hb->hbE.E[i][j], MAXHYDRO); - for (k = 0; k < MAXHYDRO; k++) - { - hb->hbE.E[i][j][k] = NULL; - } - } - } - hb->hbE.Etot = NULL; -} - -static void free_hbEmap (t_hbdata *hb) -{ - int i, j, k; - for (i = 0; i < hb->d.nrd; i++) - { - for (j = 0; j < hb->a.nra; j++) - { - for (k = 0; k < MAXHYDRO; k++) - { - sfree(hb->hbE.E[i][j][k]); - } - sfree(hb->hbE.E[i][j]); - } - sfree(hb->hbE.E[i]); - } - sfree(hb->hbE.E); - sfree(hb->hbE.Etot); -} - -static void addFramesNN(t_hbdata *hb, int frame) -{ - -#define DELTAFRAMES_HBE 10 - - int d, a, h, nframes; - - if (frame >= hb->hbE.nframes) - { - nframes = hb->hbE.nframes + DELTAFRAMES_HBE; - srenew(hb->hbE.Etot, nframes); - - for (d = 0; d < hb->d.nrd; d++) - { - for (a = 0; a < hb->a.nra; a++) - { - for (h = 0; h < hb->d.nhydro[d]; h++) - { - srenew(hb->hbE.E[d][a][h], nframes); - } - } - } - - hb->hbE.nframes += DELTAFRAMES_HBE; - } -} - -static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst, - matrix box, rvec hbox, t_donors *donors) -{ - /* d - donor atom - * a - acceptor atom - * h - hydrogen - * alpha - angle between dipoles - * x[] - atomic positions - * EEst - the type of energy estimate (see enum in hbplugin.h) - * box - the box vectors \ - * hbox - half box lengths _These two are only needed for the pbc correction - */ - - t_E E; - rvec dist; - rvec dipole[2], xmol[3], xmean[2]; - int i; - real r, realE; - - if (d == a) - { - /* Self-interaction */ - return NONSENSE_E; - } - - switch (EEst) - { - case NN_BINARY: - /* This is a simple binary existence function that sets E=1 whenever - * the distance between the oxygens is equal too or less than 0.35 nm. - */ - rvec_sub(x[d], x[a], dist); - pbc_correct_gem(dist, box, hbox); - if (norm(dist) <= 0.35) - { - E = 1; - } - else - { - E = 0; - } - break; - - case NN_1_over_r3: - /* Negative potential energy of a dipole. - * E = -cos(alpha) * 1/r^3 */ - - copy_rvec(x[d], xmol[0]); /* donor */ - copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ - - svmul(15.9994*(1/1.008), xmol[0], xmean[0]); - rvec_inc(xmean[0], xmol[1]); - rvec_inc(xmean[0], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - /* Assumes that all acceptors are also donors. */ - copy_rvec(x[a], xmol[0]); /* acceptor */ - copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */ - - - svmul(15.9994*(1/1.008), xmol[0], xmean[1]); - rvec_inc(xmean[1], xmol[1]); - rvec_inc(xmean[1], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - rvec_sub(xmean[0], xmean[1], dist); - pbc_correct_gem(dist, box, hbox); - r = norm(dist); - - realE = pow(r, -3.0); - E = (t_E)(SCALEFACTOR_E * realE); - break; - - case NN_dipole: - /* Negative potential energy of a (unpolarizable) dipole. - * E = -cos(alpha) * 1/r^3 */ - clear_rvec(dipole[1]); - clear_rvec(dipole[0]); - - copy_rvec(x[d], xmol[0]); /* donor */ - copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ - - rvec_inc(dipole[0], xmol[1]); - rvec_inc(dipole[0], xmol[2]); - for (i = 0; i < 3; i++) - { - dipole[0][i] *= 0.5; - } - rvec_dec(dipole[0], xmol[0]); - - svmul(15.9994*(1/1.008), xmol[0], xmean[0]); - rvec_inc(xmean[0], xmol[1]); - rvec_inc(xmean[0], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - /* Assumes that all acceptors are also donors. */ - copy_rvec(x[a], xmol[0]); /* acceptor */ - copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ - copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */ - - - rvec_inc(dipole[1], xmol[1]); - rvec_inc(dipole[1], xmol[2]); - for (i = 0; i < 3; i++) - { - dipole[1][i] *= 0.5; - } - rvec_dec(dipole[1], xmol[0]); - - svmul(15.9994*(1/1.008), xmol[0], xmean[1]); - rvec_inc(xmean[1], xmol[1]); - rvec_inc(xmean[1], xmol[2]); - for (i = 0; i < 3; i++) - { - xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; - } - - rvec_sub(xmean[0], xmean[1], dist); - pbc_correct_gem(dist, box, hbox); - r = norm(dist); - - double cosalpha = cos_angle(dipole[0], dipole[1]); - realE = cosalpha * pow(r, -3.0); - E = (t_E)(SCALEFACTOR_E * realE); - break; - - default: - printf("Can't do that type of energy estimate: %i\n.", EEst); - E = NONSENSE_E; - } - - return E; -} - -static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame) -{ - /* hb - hbond data structure - d - donor - a - acceptor - h - hydrogen - E - estimate of the energy - frame - the current frame. - */ - - /* Store the estimated energy */ - if (E == NONSENSE_E) - { - E = 0; - } - - hb->hbE.E[d][a][h][frame] = E; - -#pragma omp critical - { - hb->hbE.Etot[frame] += E; - } -} -#endif /* HAVE_NN_LOOPS */ - - -/* Finds -v[] in the periodicity index */ -static int findMirror(PSTYPE p, ivec v[], PSTYPE nper) -{ - PSTYPE i; - ivec u; - for (i = 0; i < nper; i++) - { - if (v[i][XX] == -(v[p][XX]) && - v[i][YY] == -(v[p][YY]) && - v[i][ZZ] == -(v[p][ZZ])) - { - return (int)i; - } - } - printf("Couldn't find mirror of [%i, %i, %i], index \n", - v[p][XX], - v[p][YY], - v[p][ZZ]); - return -1; -} - - static void add_frames(t_hbdata *hb, int nframes) { int i, j, k, l; @@ -727,64 +280,13 @@ static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb) _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE); } -static void addPshift(t_pShift *pHist, PSTYPE p, int frame) -{ - if (pHist->len == 0) - { - snew(pHist->frame, 1); - snew(pHist->p, 1); - pHist->len = 1; - pHist->frame[0] = frame; - pHist->p[0] = p; - return; - } - else - if (pHist->p[pHist->len-1] != p) - { - pHist->len++; - srenew(pHist->frame, pHist->len); - srenew(pHist->p, pHist->len); - pHist->frame[pHist->len-1] = frame; - pHist->p[pHist->len-1] = p; - } /* Otherwise, there is no transition. */ - return; -} - -static PSTYPE getPshift(t_pShift pHist, int frame) -{ - int f, i; - - if (pHist.len == 0 - || (pHist.len > 0 && pHist.frame[0] > frame)) - { - return -1; - } - - for (i = 0; i < pHist.len; i++) - { - f = pHist.frame[i]; - if (f == frame) - { - return pHist.p[i]; - } - if (f > frame) - { - return pHist.p[i-1]; - } - } - - /* It seems that frame is after the last periodic transition. Return the last periodicity. */ - return pHist.p[pHist.len-1]; -} - -static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPE p) +static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb) { int i, j, n; t_hbond *hb = hbd->hbmap[id][ia]; int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id]); int wlen = hbd->wordlen; int delta = 32*wlen; - gmx_bool bGem = hbd->bGem; if (!hb->h[0]) { @@ -822,18 +324,6 @@ static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PST if (frame >= 0) { set_hb(hbd, id, h, ia, frame, ihb); - if (bGem) - { - if (p >= hbd->per->nper) - { - gmx_fatal(FARGS, "invalid shift: p=%u, nper=%u", p, hbd->per->nper); - } - else - { - addPshift(&(hbd->per->pHist[id][ia]), p, frame); - } - - } } } @@ -917,7 +407,7 @@ static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd) static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, - int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPE p) + int frame, gmx_bool bMerge, int ihb, gmx_bool bContact) { int k, id, ia, hh; gmx_bool daSwap = FALSE; @@ -1011,7 +501,7 @@ static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, snew(hb->hbmap[id][ia]->h, hb->maxhydro); snew(hb->hbmap[id][ia]->g, hb->maxhydro); } - add_ff(hb, id, k, ia, frame, ihb, p); + add_ff(hb, id, k, ia, frame, ihb); } } @@ -1503,52 +993,30 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, rvec_sub(x[ad[i]], xshell, dshell); if (bBox) { - if (FALSE && !hb->bGem) + gmx_bool bDone = FALSE; + while (!bDone) { + bDone = TRUE; for (m = DIM-1; m >= 0 && bInShell; m--) { if (dshell[m] < -hbox[m]) { + bDone = FALSE; rvec_inc(dshell, box[m]); } - else if (dshell[m] >= hbox[m]) + if (dshell[m] >= hbox[m]) { + bDone = FALSE; dshell[m] -= 2*hbox[m]; } - /* if we're outside the cube, we're outside the sphere also! */ - if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) - { - bInShell = FALSE; - } } } - else + for (m = DIM-1; m >= 0 && bInShell; m--) { - gmx_bool bDone = FALSE; - while (!bDone) - { - bDone = TRUE; - for (m = DIM-1; m >= 0 && bInShell; m--) - { - if (dshell[m] < -hbox[m]) - { - bDone = FALSE; - rvec_inc(dshell, box[m]); - } - if (dshell[m] >= hbox[m]) - { - bDone = FALSE; - dshell[m] -= 2*hbox[m]; - } - } - } - for (m = DIM-1; m >= 0 && bInShell; m--) + /* if we're outside the cube, we're outside the sphere also! */ + if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) { - /* if we're outside the cube, we're outside the sphere also! */ - if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) - { - bInShell = FALSE; - } + bInShell = FALSE; } } } @@ -1563,10 +1031,6 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, { if (bBox) { - if (hb->bGem) - { - copy_rvec(x[ad[i]], xtemp); - } pbc_in_gridbox(x[ad[i]], box); for (m = DIM-1; m >= 0; m--) @@ -1574,10 +1038,6 @@ static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, grididx[m] = x[ad[i]][m]*invdelta[m]; grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m]; } - if (hb->bGem) - { - copy_rvec(xtemp, x[ad[i]]); /* copy back */ - } } gx = grididx[XX]; @@ -1774,7 +1234,7 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, real rcut, real r2cut, real ccut, rvec x[], gmx_bool bBox, matrix box, rvec hbox, real *d_ha, real *ang, gmx_bool bDA, int *hhh, - gmx_bool bContact, gmx_bool bMerge, PSTYPE *p) + gmx_bool bContact, gmx_bool bMerge) { int h, hh, id, ja, ihb; rvec r_da, r_ha, r_dh, r = {0, 0, 0}; @@ -1811,15 +1271,7 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, { /* return hbNo; */ daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */ } - if (hb->bGem) - { - copy_rvec(r_da, r); /* Save this for later */ - pbc_correct_gem(r_da, box, hbox); - } - else - { - pbc_correct_gem(r_da, box, hbox); - } + pbc_correct_gem(r_da, box, hbox); } rda2 = iprod(r_da, r_da); @@ -1831,11 +1283,6 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, } if (rda2 <= rc2) { - if (hb->bGem) - { - calcBoxDistance(hb->per->P, r, ri); - *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */ - } return hbHB; } else if (rda2 < r2c2) @@ -1868,12 +1315,6 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, rha2 = iprod(r_ha, r_ha); } - if (hb->bGem) - { - calcBoxDistance(hb->per->P, r, ri); - *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */ - } - if (bDA || (!bDA && (rha2 <= rc2))) { rvec_sub(x[d], x[hh], r_dh); @@ -1907,28 +1348,17 @@ static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, } } -/* Fixed previously undiscovered bug in the merge - code, where the last frame of each hbond disappears. - - Erik Marklund, June 1, 2006 */ -/* Added the following arguments: - * ptmp[] - temporary periodicity hisory - * a1 - identity of first acceptor/donor - * a2 - identity of second acceptor/donor - * - Erik Marklund, FEB 20 2010 */ - /* Merging is now done on the fly, so do_merge is most likely obsolete now. * Will do some more testing before removing the function entirely. * - Erik Marklund, MAY 10 2010 */ static void do_merge(t_hbdata *hb, int ntmp, - unsigned int htmp[], unsigned int gtmp[], PSTYPE ptmp[], - t_hbond *hb0, t_hbond *hb1, int a1, int a2) + unsigned int htmp[], unsigned int gtmp[], + t_hbond *hb0, t_hbond *hb1) { /* Here we need to make sure we're treating periodicity in * the right way for the geminate recombination kinetics. */ int m, mm, n00, n01, nn0, nnframes; - PSTYPE pm; - t_pShift *pShift; /* Decide where to start from when merging */ n00 = hb0->n0; @@ -1940,7 +1370,6 @@ static void do_merge(t_hbdata *hb, int ntmp, { htmp[m] = 0; gtmp[m] = 0; - ptmp[m] = 0; } /* Fill tmp arrays with values due to first HB */ /* Once again '<' had to be replaced with '<=' @@ -1951,28 +1380,11 @@ static void do_merge(t_hbdata *hb, int ntmp, { mm = m+n00-nn0; htmp[mm] = is_hb(hb0->h[0], m); - if (hb->bGem) - { - pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0); - if (pm > hb->per->nper) - { - gmx_fatal(FARGS, "Illegal shift!"); - } - else - { - ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/ - } - } } - /* If we're doing geminate recompbination we usually don't need the distances. - * Let's save some memory and time. */ - if (TRUE || !hb->bGem || hb->per->gemtype == gemAD) + for (m = 0; (m <= hb0->nframes); m++) { - for (m = 0; (m <= hb0->nframes); m++) - { - mm = m+n00-nn0; - gtmp[mm] = is_hb(hb0->g[0], m); - } + mm = m+n00-nn0; + gtmp[mm] = is_hb(hb0->g[0], m); } /* Next HB */ for (m = 0; (m <= hb1->nframes); m++) @@ -1980,20 +1392,6 @@ static void do_merge(t_hbdata *hb, int ntmp, mm = m+n01-nn0; htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m); gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m); - if (hb->bGem /* && ptmp[mm] != 0 */) - { - - /* If this hbond has been seen before with donor and acceptor swapped, - * then we need to find the mirrored (*-1) periodicity vector to truely - * merge the hbond history. */ - pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper); - /* Store index of mirror */ - if (pm > hb->per->nper) - { - gmx_fatal(FARGS, "Illegal shift!"); - } - ptmp[mm] = pm; - } } /* Reallocate target array */ if (nnframes > hb0->maxframes) @@ -2001,20 +1399,12 @@ static void do_merge(t_hbdata *hb, int ntmp, srenew(hb0->h[0], 4+nnframes/hb->wordlen); srenew(hb0->g[0], 4+nnframes/hb->wordlen); } - if (NULL != hb->per->pHist) - { - clearPshift(&(hb->per->pHist[a1][a2])); - } /* Copy temp array to target array */ for (m = 0; (m <= nnframes); m++) { _set_hb(hb0->h[0], m, htmp[m]); _set_hb(hb0->g[0], m, gtmp[m]); - if (hb->bGem) - { - addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0); - } } /* Set scalar variables */ @@ -2022,14 +1412,10 @@ static void do_merge(t_hbdata *hb, int ntmp, hb0->maxframes = nnframes; } -/* Added argument bContact for nicer output. - * Erik Marklund, June 29, 2006 - */ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) { int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp; unsigned int *htmp, *gtmp; - PSTYPE *ptmp; t_hbond *hb0, *hb1; inrnew = hb->nrhb; @@ -2041,7 +1427,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) ntmp = 2*hb->max_frames; snew(gtmp, ntmp); snew(htmp, ntmp); - snew(ptmp, ntmp); for (i = 0; (i < hb->d.nrd); i++) { fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd); @@ -2058,7 +1443,7 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) hb1 = hb->hbmap[jj][ii]; if (hb0 && hb1 && ISHB(hb0->history[0]) && ISHB(hb1->history[0])) { - do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j); + do_merge(hb, ntmp, htmp, gtmp, hb0, hb1); if (ISHB(hb1->history[0])) { inrnew--; @@ -2078,10 +1463,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) } sfree(hb1->h[0]); sfree(hb1->g[0]); - if (hb->bGem) - { - clearPshift(&(hb->per->pHist[jj][ii])); - } hb1->h[0] = NULL; hb1->g[0] = NULL; hb1->history[0] = hbNo; @@ -2096,7 +1477,6 @@ static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) hb->nrdist = indnew; sfree(gtmp); sfree(htmp); - sfree(ptmp); } static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) @@ -2128,12 +1508,6 @@ static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) fprintf(fp, " %8d\n", nbtot); } -/* Added argument bContact in do_hblife(...). Also - * added support for -contact in function body. - * - Erik Marklund, May 31, 2006 */ -/* Changed the contact code slightly. - * - Erik Marklund, June 29, 2006 - */ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact, const output_env_t oenv) { @@ -2185,10 +1559,6 @@ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bC ohb = 0; j0 = 0; - /* Changed '<' into '<=' below, just like I - did in the hbm-output-loop in the main code. - - Erik Marklund, May 31, 2006 - */ for (j = 0; (j <= hbh->nframes); j++) { ihb = is_hb(h[nh], j); @@ -2256,9 +1626,6 @@ static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bC sfree(histo); } -/* Changed argument bMerge into oneHB to handle contacts properly. - * - Erik Marklund, June 29, 2006 - */ static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump) { FILE *fp; @@ -2549,17 +1916,10 @@ static void normalizeACF(real *ct, real *gt, int nhb, int len) } } -/* Added argument bContact in do_hbac(...). Also - * added support for -contact in the actual code. - * - Erik Marklund, May 31, 2006 */ -/* Changed contact code and added argument R2 - * - Erik Marklund, June 29, 2006 - */ static void do_hbac(const char *fn, t_hbdata *hb, int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start, real temp, gmx_bool R2, const output_env_t oenv, - const char *gemType, int nThreads, - const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit) + int nThreads) { FILE *fp; int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets; @@ -2585,13 +1945,9 @@ static void do_hbac(const char *fn, t_hbdata *hb, unsigned int **h = NULL, **g = NULL; int nh, nhbonds, nhydro, ngh; t_hbond *hbh; - PSTYPE p, *pfound = NULL, np; - t_pShift *pHist; int *ptimes = NULL, *poff = NULL, anhb, n0, mMax = INT_MIN; - real **rHbExGem = NULL; gmx_bool c; int acType; - t_E *E; double *ctdouble, *timedouble, *fittedct; double fittolerance = 0.1; int *dondata = NULL, thisThread; @@ -2608,44 +1964,8 @@ static void do_hbac(const char *fn, t_hbdata *hb, printf("Doing autocorrelation "); - /* Decide what kind of ACF calculations to do. */ - if (NN > NN_NONE && NN < NN_NR) - { -#ifdef HAVE_NN_LOOPS - acType = AC_NN; - printf("using the energy estimate.\n"); -#else - acType = AC_NONE; - printf("Can't do the NN-loop. Yet.\n"); -#endif - } - else if (hb->bGem) - { - acType = AC_GEM; - printf("according to the reversible geminate recombination model by Omer Markowitch.\n"); - - nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0); - snew(legGem, nSets); - for (i = 0; i < nSets; i++) - { - snew(legGem[i], 128); - } - sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType); - if (bBallistic) - { - sprintf(legGem[1], "Ac'(t)"); - } - if (bGemFit) - { - sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType); - } - - } - else - { - acType = AC_LUZAR; - printf("according to the theory of Luzar and Chandler.\n"); - } + acType = AC_LUZAR; + printf("according to the theory of Luzar and Chandler.\n"); fflush(stdout); /* build hbexist matrix in reals for autocorr */ @@ -2698,581 +2018,173 @@ static void do_hbac(const char *fn, t_hbdata *hb, } - /* Build the ACF according to acType */ - switch (acType) - { - - case AC_NN: -#ifdef HAVE_NN_LOOPS - /* Here we're using the estimated energy for the hydrogen bonds. */ - snew(ct, nn); - -#pragma omp parallel \ - private(i, j, k, nh, E, rhbex, thisThread) \ - default(shared) - { -#pragma omp barrier - thisThread = gmx_omp_get_thread_num(); - rhbex = NULL; - - snew(rhbex, n2); - memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */ + /* Build the ACF */ + snew(rhbex, 2*n2); + snew(ct, 2*n2); + snew(gt, 2*n2); + snew(ht, 2*n2); + snew(ght, 2*n2); + snew(dght, 2*n2); -#pragma omp barrier -#pragma omp for schedule (dynamic) - for (i = 0; i < hb->d.nrd; i++) /* loop over donors */ - { - if (bOMP) - { -#pragma omp critical - { - dondata[thisThread] = i; - parallel_print(dondata, nThreads); - } - } - else - { - fprintf(stderr, "\r %i", i); - } + snew(kt, nn); + snew(cct, nn); - for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */ - { - for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */ - { - E = hb->hbE.E[i][j][nh]; - if (E != NULL) - { - for (k = 0; k < nframes; k++) - { - if (E[k] != NONSENSE_E) - { - rhbex[k] = (real)E[k]; - } - } - - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0, 1); -#pragma omp critical - { - for (k = 0; (k < nn); k++) - { - ct[k] += rhbex[k]; - } - } - } - } /* k loop */ - } /* j loop */ - } /* i loop */ - sfree(rhbex); -#pragma omp barrier - } - - if (bOMP) - { - sfree(dondata); - } - normalizeACF(ct, NULL, 0, nn); - snew(ctdouble, nn); - snew(timedouble, nn); - for (j = 0; j < nn; j++) - { - timedouble[j] = (double)(hb->time[j]); - ctdouble[j] = (double)(ct[j]); - } - - /* Remove ballistic term */ - /* Ballistic component removal and fitting to the reversible geminate recombination model - * will be taken out for the time being. First of all, one can remove the ballistic - * component with g_analyze afterwards. Secondly, and more importantly, there are still - * problems with the robustness of the fitting to the model. More work is needed. - * A third reason is that we're currently using gsl for this and wish to reduce dependence - * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with - * a BSD-licence that can do the job. - * - * - Erik Marklund, June 18 2010. - */ -/* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ -/* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ -/* else */ -/* printf("\nNumber of data points is less than the number of parameters to fit\n." */ -/* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ - - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)"); - xvgr_legend(fp, asize(legNN), legNN); + for (i = 0; (i < hb->d.nrd); i++) + { + for (k = 0; (k < hb->a.nra); k++) + { + nhydro = 0; + hbh = hb->hbmap[i][k]; - for (j = 0; (j < nn); j++) + if (hbh) { - fprintf(fp, "%10g %10g %10g\n", - hb->time[j]-hb->time[0], - ct[j], - ctdouble[j]); - } - xvgrclose(fp); - sfree(ct); - sfree(ctdouble); - sfree(timedouble); -#endif /* HAVE_NN_LOOPS */ - break; /* case AC_NN */ - - case AC_GEM: - snew(ct, 2*n2); - memset(ct, 0, 2*n2*sizeof(real)); -#ifndef GMX_OPENMP - fprintf(stderr, "Donor:\n"); -#define __ACDATA ct -#else -#define __ACDATA p_ct -#endif - -#pragma omp parallel \ - private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \ - pfound, poff, rHbExGem, p, ihb, mMax, \ - thisThread, p_ct) \ - default(shared) - { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */ - h = NULL; - g = NULL; - thisThread = gmx_omp_get_thread_num(); - snew(h, hb->maxhydro); - snew(g, hb->maxhydro); - mMax = INT_MIN; - rHbExGem = NULL; - poff = NULL; - pfound = NULL; - p_ct = NULL; - snew(p_ct, 2*n2); - memset(p_ct, 0, 2*n2*sizeof(real)); - - /* I'm using a chunk size of 1, since I expect \ - * the overhead to be really small compared \ - * to the actual calculations \ */ -#pragma omp for schedule(dynamic,1) nowait - for (i = 0; i < hb->d.nrd; i++) + if (bMerge || bContact) { - - if (bOMP) + if (ISHB(hbh->history[0])) { -#pragma omp critical - { - dondata[thisThread] = i; - parallel_print(dondata, nThreads); - } - } - else - { - fprintf(stderr, "\r %i", i); + h[0] = hbh->h[0]; + g[0] = hbh->g[0]; + nhydro = 1; } - for (k = 0; k < hb->a.nra; k++) - { - for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) - { - hbh = hb->hbmap[i][k]; - if (hbh) - { - /* Note that if hb->per->gemtype==gemDD, then distances will be stored in - * hb->hbmap[d][a].h array anyway, because the contact flag will be set. - * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */ - pHist = &(hb->per->pHist[i][k]); - if (ISHB(hbh->history[nh]) && pHist->len != 0) - { - - { - h[nh] = hbh->h[nh]; - g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL; - } - n0 = hbh->n0; - nf = hbh->nframes; - /* count the number of periodic shifts encountered and store - * them in separate arrays. */ - np = 0; - for (j = 0; j < pHist->len; j++) - { - p = pHist->p[j]; - for (m = 0; m <= np; m++) - { - if (m == np) /* p not recognized in list. Add it and set up new array. */ - { - np++; - if (np > hb->per->nper) - { - gmx_fatal(FARGS, "Too many pshifts. Something's utterly wrong here."); - } - if (m >= mMax) /* Extend the arrays. - * Doing it like this, using mMax to keep track of the sizes, - * eleviates the need for freeing and re-allocating the arrays - * when taking on the next donor-acceptor pair */ - { - mMax = m; - srenew(pfound, np); /* The list of found periodic shifts. */ - srenew(rHbExGem, np); /* The hb existence functions (-aver_hb). */ - snew(rHbExGem[m], 2*n2); - srenew(poff, np); - } - - { - if (rHbExGem != NULL && rHbExGem[m] != NULL) - { - /* This must be done, as this array was most likey - * used to store stuff in some previous iteration. */ - memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2)); - } - else - { - fprintf(stderr, "rHbExGem not initialized! m = %i\n", m); - } - } - pfound[m] = p; - poff[m] = -1; - - break; - } /* m==np */ - if (p == pfound[m]) - { - break; - } - } /* m: Loop over found shifts */ - } /* j: Loop over shifts */ - - /* Now unpack and disentangle the existence funtions. */ - for (j = 0; j < nf; j++) - { - /* i: donor, - * k: acceptor - * nh: hydrogen - * j: time - * p: periodic shift - * pfound: list of periodic shifts found for this pair. - * poff: list of frame offsets; that is, the first - * frame a hbond has a particular periodic shift. */ - p = getPshift(*pHist, j+n0); - if (p != -1) - { - for (m = 0; m < np; m++) - { - if (pfound[m] == p) - { - break; - } - if (m == (np-1)) - { - gmx_fatal(FARGS, "Shift not found, but must be there."); - } - } - - ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE : is_hb(g[nh], j)); - if (ihb) - { - if (poff[m] == -1) - { - poff[m] = j; /* Here's where the first hbond with shift p is, - * relative to the start of h[0].*/ - } - if (j < poff[m]) - { - gmx_fatal(FARGS, "j 0 && n0+poff[m] < nn /* && m==0 */) - { - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); - for (j = 0; (j < nn); j++) - { - __ACDATA[j] += rHbExGem[m][j]; - } - } - } /* Building of ac. */ - } /* if (ISHB(...*/ - } /* if (hbh) */ - } /* hydrogen loop */ - } /* acceptor loop */ - } /* donor loop */ - - for (m = 0; m <= mMax; m++) - { - sfree(rHbExGem[m]); } - sfree(pfound); - sfree(poff); - sfree(rHbExGem); - - sfree(h); - sfree(g); - - if (bOMP) + else { -#pragma omp critical + for (m = 0; (m < hb->maxhydro); m++) { - for (i = 0; i < nn; i++) + if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) { - ct[i] += p_ct[i]; + g[nhydro] = hbh->g[m]; + h[nhydro] = hbh->h[m]; + nhydro++; } } - sfree(p_ct); } - } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */ - if (bOMP) - { - sfree(dondata); - } - - normalizeACF(ct, NULL, 0, nn); - - fprintf(stderr, "\n\nACF successfully calculated.\n"); - - /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */ - - snew(ctdouble, nn); - snew(timedouble, nn); - snew(fittedct, nn); - - for (j = 0; j < nn; j++) - { - timedouble[j] = (double)(hb->time[j]); - ctdouble[j] = (double)(ct[j]); - } - - /* Remove ballistic term */ - /* Ballistic component removal and fitting to the reversible geminate recombination model - * will be taken out for the time being. First of all, one can remove the ballistic - * component with g_analyze afterwards. Secondly, and more importantly, there are still - * problems with the robustness of the fitting to the model. More work is needed. - * A third reason is that we're currently using gsl for this and wish to reduce dependence - * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with - * a BSD-licence that can do the job. - * - * - Erik Marklund, June 18 2010. - */ -/* if (bBallistic) { */ -/* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ -/* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ -/* else */ -/* printf("\nNumber of data points is less than the number of parameters to fit\n." */ -/* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ -/* } */ -/* if (bGemFit) */ -/* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */ - - - if (bContact) - { - fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - else - { - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - xvgr_legend(fp, asize(legGem), (const char**)legGem, oenv); - - for (j = 0; (j < nn); j++) - { - fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]); - if (bBallistic) - { - fprintf(fp, " %10g", ctdouble[j]); - } - if (bGemFit) - { - fprintf(fp, " %10g", fittedct[j]); - } - fprintf(fp, "\n"); - } - xvgrclose(fp); - - sfree(ctdouble); - sfree(timedouble); - sfree(fittedct); - sfree(ct); - - break; /* case AC_GEM */ - - case AC_LUZAR: - snew(rhbex, 2*n2); - snew(ct, 2*n2); - snew(gt, 2*n2); - snew(ht, 2*n2); - snew(ght, 2*n2); - snew(dght, 2*n2); - - snew(kt, nn); - snew(cct, nn); - - for (i = 0; (i < hb->d.nrd); i++) - { - for (k = 0; (k < hb->a.nra); k++) + nf = hbh->nframes; + for (nh = 0; (nh < nhydro); nh++) { - nhydro = 0; - hbh = hb->hbmap[i][k]; - - if (hbh) + int nrint = bContact ? hb->nrdist : hb->nrhb; + if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) + { + fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint); + } + nhbonds++; + for (j = 0; (j < nframes); j++) { - if (bMerge || bContact) + if (j <= nf) { - if (ISHB(hbh->history[0])) - { - h[0] = hbh->h[0]; - g[0] = hbh->g[0]; - nhydro = 1; - } + ihb = is_hb(h[nh], j); + idist = is_hb(g[nh], j); } else { - for (m = 0; (m < hb->maxhydro); m++) - { - if (bContact ? ISDIST(hbh->history[m]) : ISHB(hbh->history[m])) - { - g[nhydro] = hbh->g[m]; - h[nhydro] = hbh->h[m]; - nhydro++; - } - } + ihb = idist = 0; } - - nf = hbh->nframes; - for (nh = 0; (nh < nhydro); nh++) + rhbex[j] = ihb; + /* For contacts: if a second cut-off is provided, use it, + * otherwise use g(t) = 1-h(t) */ + if (!R2 && bContact) { - int nrint = bContact ? hb->nrdist : hb->nrhb; - if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) - { - fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint); - } - nhbonds++; - for (j = 0; (j < nframes); j++) - { - /* Changed '<' into '<=' below, just like I did in - the hbm-output-loop in the gmx_hbond() block. - - Erik Marklund, May 31, 2006 */ - if (j <= nf) - { - ihb = is_hb(h[nh], j); - idist = is_hb(g[nh], j); - } - else - { - ihb = idist = 0; - } - rhbex[j] = ihb; - /* For contacts: if a second cut-off is provided, use it, - * otherwise use g(t) = 1-h(t) */ - if (!R2 && bContact) - { - gt[j] = 1-ihb; - } - else - { - gt[j] = idist*(1-ihb); - } - ht[j] = rhbex[j]; - nhb += ihb; - } - + gt[j] = 1-ihb; + } + else + { + gt[j] = idist*(1-ihb); + } + ht[j] = rhbex[j]; + nhb += ihb; + } - /* The autocorrelation function is normalized after summation only */ - low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], - eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); + /* The autocorrelation function is normalized after summation only */ + low_do_autocorr(NULL, oenv, NULL, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], + eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0); - /* Cross correlation analysis for thermodynamics */ - for (j = nframes; (j < n2); j++) - { - ht[j] = 0; - gt[j] = 0; - } + /* Cross correlation analysis for thermodynamics */ + for (j = nframes; (j < n2); j++) + { + ht[j] = 0; + gt[j] = 0; + } - cross_corr(n2, ht, gt, dght); + cross_corr(n2, ht, gt, dght); - for (j = 0; (j < nn); j++) - { - ct[j] += rhbex[j]; - ght[j] += dght[j]; - } - } + for (j = 0; (j < nn); j++) + { + ct[j] += rhbex[j]; + ght[j] += dght[j]; } } } - fprintf(stderr, "\n"); - sfree(h); - sfree(g); - normalizeACF(ct, ght, nhb, nn); - - /* Determine tail value for statistics */ - tail = 0; - tail2 = 0; - for (j = nn/2; (j < nn); j++) - { - tail += ct[j]; - tail2 += ct[j]*ct[j]; - } - tail /= (nn - nn/2); - tail2 /= (nn - nn/2); - dtail = sqrt(tail2-tail*tail); + } + } + fprintf(stderr, "\n"); + sfree(h); + sfree(g); + normalizeACF(ct, ght, nhb, nn); - /* Check whether the ACF is long enough */ - if (dtail > tol) - { - printf("\nWARNING: Correlation function is probably not long enough\n" - "because the standard deviation in the tail of C(t) > %g\n" - "Tail value (average C(t) over second half of acf): %g +/- %g\n", - tol, tail, dtail); - } - for (j = 0; (j < nn); j++) - { - cct[j] = ct[j]; - ct[j] = (cct[j]-tail)/(1-tail); - } - /* Compute negative derivative k(t) = -dc(t)/dt */ - compute_derivative(nn, hb->time, ct, kt); - for (j = 0; (j < nn); j++) - { - kt[j] = -kt[j]; - } + /* Determine tail value for statistics */ + tail = 0; + tail2 = 0; + for (j = nn/2; (j < nn); j++) + { + tail += ct[j]; + tail2 += ct[j]*ct[j]; + } + tail /= (nn - nn/2); + tail2 /= (nn - nn/2); + dtail = sqrt(tail2-tail*tail); + /* Check whether the ACF is long enough */ + if (dtail > tol) + { + printf("\nWARNING: Correlation function is probably not long enough\n" + "because the standard deviation in the tail of C(t) > %g\n" + "Tail value (average C(t) over second half of acf): %g +/- %g\n", + tol, tail, dtail); + } + for (j = 0; (j < nn); j++) + { + cct[j] = ct[j]; + ct[j] = (cct[j]-tail)/(1-tail); + } + /* Compute negative derivative k(t) = -dc(t)/dt */ + compute_derivative(nn, hb->time, ct, kt); + for (j = 0; (j < nn); j++) + { + kt[j] = -kt[j]; + } - if (bContact) - { - fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - else - { - fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); - } - xvgr_legend(fp, asize(legLuzar), legLuzar, oenv); + if (bContact) + { + fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); + } + else + { + fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); + } + xvgr_legend(fp, asize(legLuzar), legLuzar, oenv); - for (j = 0; (j < nn); j++) - { - fprintf(fp, "%10g %10g %10g %10g %10g\n", - hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); - } - xvgrclose(fp); - - analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL, - fit_start, temp); - - do_view(oenv, fn, NULL); - sfree(rhbex); - sfree(ct); - sfree(gt); - sfree(ht); - sfree(ght); - sfree(dght); - sfree(cct); - sfree(kt); - /* sfree(h); */ -/* sfree(g); */ - - break; /* case AC_LUZAR */ - - default: - gmx_fatal(FARGS, "Unrecognized type of ACF-calulation. acType = %i.", acType); - } /* switch (acType) */ + + for (j = 0; (j < nn); j++) + { + fprintf(fp, "%10g %10g %10g %10g %10g\n", + hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); + } + xvgrclose(fp); + + analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL, + fit_start, temp); + + do_view(oenv, fn, NULL); + sfree(rhbex); + sfree(ct); + sfree(gt); + sfree(ht); + sfree(ght); + sfree(dght); + sfree(cct); + sfree(kt); } static void init_hbframe(t_hbdata *hb, int nframes, real t) @@ -3286,24 +2198,6 @@ static void init_hbframe(t_hbdata *hb, int nframes, real t) { hb->nhx[nframes][i] = 0; } - /* Loop invalidated */ - if (hb->bHBmap && 0) - { - for (i = 0; (i < hb->d.nrd); i++) - { - for (j = 0; (j < hb->a.nra); j++) - { - for (m = 0; (m < hb->maxhydro); m++) - { - if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m]) - { - set_hb(hb, i, m, j, nframes, HB_NO); - } - } - } - } - } - /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/ } static FILE *open_donor_properties_file(const char *fn, @@ -3383,10 +2277,7 @@ static void dump_hbmap(t_hbdata *hb, fprintf(fp, " %4d", index[grp][i]+1); } fprintf(fp, "\n"); - /* - Added -contact support below. - - Erik Marklund, May 29, 2006 - */ + if (!bContact) { fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]); @@ -3579,7 +2470,7 @@ int gmx_hbond(int argc, char *argv[]) static int nDump = 0, nFitPoints = 100; static int nThreads = 0, nBalExp = 4; - static gmx_bool bContact = FALSE, bBallistic = FALSE, bGemFit = FALSE; + static gmx_bool bContact = FALSE; static real logAfterTime = 10, gemBallistic = 0.2; /* ps */ static const char *NNtype[] = {NULL, "none", "binary", "oneOverR3", "dipole", NULL}; @@ -3616,10 +2507,6 @@ int gmx_hbond(int argc, char *argv[]) "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" }, { "-merge", FALSE, etBOOL, {&bMerge}, "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }, - { "-geminate", FALSE, etENUM, {gemType}, - "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."}, - { "-diff", FALSE, etREAL, {&D}, - "HIDDENDffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."}, #ifdef GMX_OPENMP { "-nthreads", FALSE, etINT, {&nThreads}, "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"}, @@ -3681,13 +2568,10 @@ int gmx_hbond(int argc, char *argv[]) ivec ngrid; unsigned char *datable; output_env_t oenv; - int gemmode, NN; - PSTYPE peri = 0; - t_E E; + int NN; int ii, jj, hh, actual_nThreads; int threadNr = 0; - gmx_bool bGem, bNN, bParallel; - t_gemParams *params = NULL; + gmx_bool bParallel; gmx_bool bEdge_yjj, bEdge_xjj, bOMP; t_hbdata **p_hb = NULL; /* one per thread, then merge after the frame loop */ @@ -3708,79 +2592,6 @@ int gmx_hbond(int argc, char *argv[]) return 0; } - /* NN-loop? If so, what estimator to use ?*/ - NN = 1; - /* Outcommented for now DvdS 2010-07-13 - while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0) - NN++; - if (NN == NN_NR) - gmx_fatal(FARGS, "Invalid NN-loop type."); - */ - bNN = FALSE; - for (i = 2; bNN == FALSE && i < NN_NR; i++) - { - bNN = bNN || NN == i; - } - - if (NN > NN_NONE && bMerge) - { - bMerge = FALSE; - } - - /* geminate recombination? If so, which flavor? */ - gemmode = 1; - while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0) - { - gemmode++; - } - if (gemmode == gemNR) - { - gmx_fatal(FARGS, "Invalid recombination type."); - } - - bGem = FALSE; - for (i = 2; bGem == FALSE && i < gemNR; i++) - { - bGem = bGem || gemmode == i; - } - - if (bGem) - { - printf("Geminate recombination: %s\n", gemType[gemmode]); - if (bContact) - { - if (gemmode != gemDD) - { - printf("Turning off -contact option...\n"); - bContact = FALSE; - } - } - else - { - if (gemmode == gemDD) - { - printf("Turning on -contact option...\n"); - bContact = TRUE; - } - } - if (bMerge) - { - if (gemmode == gemAA) - { - printf("Turning off -merge option...\n"); - bMerge = FALSE; - } - } - else - { - if (gemmode != gemAA) - { - printf("Turning on -merge option...\n"); - bMerge = TRUE; - } - } - } - /* process input */ bSelected = FALSE; ccut = cos(acut*DEG2RAD); @@ -3801,8 +2612,7 @@ int gmx_hbond(int argc, char *argv[]) bHBmap = (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm) || opt2bSet("-hbn", NFILE, fnm) || - opt2bSet("-hbm", NFILE, fnm) || - bGem); + opt2bSet("-hbm", NFILE, fnm)); if (opt2bSet("-nhbdist", NFILE, fnm)) { @@ -3812,7 +2622,7 @@ int gmx_hbond(int argc, char *argv[]) xvgr_legend(fpnhb, asize(leg), leg, oenv); } - hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode); + hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact); /* get topology */ read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top); @@ -3848,7 +2658,7 @@ int gmx_hbond(int argc, char *argv[]) /* Should this be here ? */ snew(hb->d.dptr, top.atoms.nr); snew(hb->a.aptr, top.atoms.nr); - add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri); + add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact); } printf("Analyzing %d selected hydrogen bonds from '%s'\n", isize[0], grpnames[0]); @@ -3944,21 +2754,6 @@ int gmx_hbond(int argc, char *argv[]) printf("done.\n"); } -#ifdef HAVE_NN_LOOPS - if (bNN) - { - mk_hbEmap(hb, 0); - } -#endif - - if (bGem) - { - printf("Making per structure..."); - /* Generate hbond data structure */ - mk_per(hb); - printf("done.\n"); - } - /* check input */ bStop = FALSE; if (hb->d.nrd + hb->a.nra == 0) @@ -4022,11 +2817,6 @@ int gmx_hbond(int argc, char *argv[]) snew(adist, nabin+1); snew(rdist, nrbin+1); - if (bGem && !bBox) - { - gmx_fatal(FARGS, "Can't do geminate recombination without periodic box."); - } - bParallel = FALSE; #ifndef GMX_OPENMP @@ -4077,7 +2867,6 @@ int gmx_hbond(int argc, char *argv[]) p_hb[i]->bHBmap = hb->bHBmap; p_hb[i]->bDAnr = hb->bDAnr; - p_hb[i]->bGem = hb->bGem; p_hb[i]->wordlen = hb->wordlen; p_hb[i]->nframes = hb->nframes; p_hb[i]->maxhydro = hb->maxhydro; @@ -4086,11 +2875,6 @@ int gmx_hbond(int argc, char *argv[]) p_hb[i]->a = hb->a; p_hb[i]->hbmap = hb->hbmap; p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */ - p_hb[i]->per = hb->per; - -#ifdef HAVE_NN_LOOPS - p_hb[i]->hbE = hb->hbE; -#endif p_hb[i]->nrhb = 0; p_hb[i]->nrdist = 0; @@ -4102,9 +2886,9 @@ int gmx_hbond(int argc, char *argv[]) #pragma omp parallel \ firstprivate(i) \ - private(j, h, ii, jj, hh, E, \ + private(j, h, ii, jj, hh, \ xi, yi, zi, xj, yj, zj, threadNr, \ - dist, ang, peri, icell, jcell, \ + dist, ang, icell, jcell, \ grp, ogrp, ai, aj, xjj, yjj, zjj, \ xk, yk, zk, ihb, id, resdist, \ xkk, ykk, zkk, kcell, ak, k, bTric, \ @@ -4147,234 +2931,179 @@ int gmx_hbond(int argc, char *argv[]) p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */ } - if (bNN) + if (bSelected) { -#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */ - /* Loop over all atom pairs and estimate interaction energy */ #pragma omp single { - addFramesNN(hb, nframes); - } - -#pragma omp barrier -#pragma omp for schedule(dynamic) - for (i = 0; i < hb->d.nrd; i++) - { - for (j = 0; j < hb->a.nra; j++) + /* Do not parallelize this just yet. */ + /* int ii; */ + for (ii = 0; (ii < nsel); ii++) { - for (h = 0; - h < (bContact ? 1 : hb->d.nhydro[i]); - h++) - { - if (i == hb->d.nrd || j == hb->a.nra) - { - gmx_fatal(FARGS, "out of bounds"); - } - - /* Get the real atom ids */ - ii = hb->d.don[i]; - jj = hb->a.acc[j]; - hh = hb->d.hydro[i][h]; + int dd = index[0][i]; + int aa = index[0][i+2]; + /* int */ hh = index[0][i+1]; + ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, + hbox, &dist, &ang, bDA, &h, bContact, bMerge); - /* Estimate the energy from the geometry */ - E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d)); - /* Store the energy */ - storeHbEnergy(hb, i, j, h, E, nframes); + if (ihb) + { + /* add to index if not already there */ + /* Add a hbond */ + add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact); } } - } -#endif /* HAVE_NN_LOOPS */ - } /* if (bNN)*/ + } /* omp single */ + } /* if (bSelected) */ else { - if (bSelected) - { - -#pragma omp single - { - /* Do not parallelize this just yet. */ - /* int ii; */ - for (ii = 0; (ii < nsel); ii++) - { - int dd = index[0][i]; - int aa = index[0][i+2]; - /* int */ hh = index[0][i+1]; - ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, - hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); - - if (ihb) - { - /* add to index if not already there */ - /* Add a hbond */ - add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri); - } - } - } /* omp single */ - } /* if (bSelected) */ - else + /* The outer grid loop will have to do for now. */ +#pragma omp for schedule(dynamic) + for (xi = 0; xi < ngrid[XX]; xi++) { - -#pragma omp single + for (yi = 0; (yi < ngrid[YY]); yi++) { - if (bGem) + for (zi = 0; (zi < ngrid[ZZ]); zi++) { - calcBoxProjection(box, hb->per->P); - } - /* loop over all gridcells (xi,yi,zi) */ - /* Removed confusing macro, DvdS 27/12/98 */ - - } - /* The outer grid loop will have to do for now. */ -#pragma omp for schedule(dynamic) - for (xi = 0; xi < ngrid[XX]; xi++) - { - for (yi = 0; (yi < ngrid[YY]); yi++) - { - for (zi = 0; (zi < ngrid[ZZ]); zi++) + /* loop over donor groups gr0 (always) and gr1 (if necessary) */ + for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) { + icell = &(grid[zi][yi][xi].d[grp]); - /* loop over donor groups gr0 (always) and gr1 (if necessary) */ - for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) + if (bTwo) + { + ogrp = 1-grp; + } + else { - icell = &(grid[zi][yi][xi].d[grp]); + ogrp = grp; + } - if (bTwo) - { - ogrp = 1-grp; - } - else - { - ogrp = grp; - } + /* loop over all hydrogen atoms from group (grp) + * in this gridcell (icell) + */ + for (ai = 0; (ai < icell->nr); ai++) + { + i = icell->atoms[ai]; - /* loop over all hydrogen atoms from group (grp) - * in this gridcell (icell) - */ - for (ai = 0; (ai < icell->nr); ai++) + /* loop over all adjacent gridcells (xj,yj,zj) */ + for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE); + zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE); + zjj++) { - i = icell->atoms[ai]; - - /* loop over all adjacent gridcells (xj,yj,zj) */ - for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE); - zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE); - zjj++) + zj = grid_mod(zjj, ngrid[ZZ]); + bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1); + for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj); + yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj); + yjj++) { - zj = grid_mod(zjj, ngrid[ZZ]); - bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1); - for (yjj = grid_loop_begin(ngrid[YY], yi, bTric, bEdge_yjj); - yjj <= grid_loop_end(ngrid[YY], yi, bTric, bEdge_yjj); - yjj++) + yj = grid_mod(yjj, ngrid[YY]); + bEdge_xjj = + (yj == 0) || (yj == ngrid[YY] - 1) || + (zj == 0) || (zj == ngrid[ZZ] - 1); + for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj); + xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj); + xjj++) { - yj = grid_mod(yjj, ngrid[YY]); - bEdge_xjj = - (yj == 0) || (yj == ngrid[YY] - 1) || - (zj == 0) || (zj == ngrid[ZZ] - 1); - for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj); - xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj); - xjj++) + xj = grid_mod(xjj, ngrid[XX]); + jcell = &(grid[zj][yj][xj].a[ogrp]); + /* loop over acceptor atoms from other group (ogrp) + * in this adjacent gridcell (jcell) + */ + for (aj = 0; (aj < jcell->nr); aj++) { - xj = grid_mod(xjj, ngrid[XX]); - jcell = &(grid[zj][yj][xj].a[ogrp]); - /* loop over acceptor atoms from other group (ogrp) - * in this adjacent gridcell (jcell) - */ - for (aj = 0; (aj < jcell->nr); aj++) - { - j = jcell->atoms[aj]; + j = jcell->atoms[aj]; - /* check if this once was a h-bond */ - peri = -1; - ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, - hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); + /* check if this once was a h-bond */ + ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, + hbox, &dist, &ang, bDA, &h, bContact, bMerge); - if (ihb) - { - /* add to index if not already there */ - /* Add a hbond */ - add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri); + if (ihb) + { + /* add to index if not already there */ + /* Add a hbond */ + add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact); - /* make angle and distance distributions */ - if (ihb == hbHB && !bContact) + /* make angle and distance distributions */ + if (ihb == hbHB && !bContact) + { + if (dist > rcut) { - if (dist > rcut) + gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist); + } + ang *= RAD2DEG; + __ADIST[(int)( ang/abin)]++; + __RDIST[(int)(dist/rbin)]++; + if (!bTwo) + { + int id, ia; + if ((id = donor_index(&hb->d, grp, i)) == NOTSET) { - gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist); + gmx_fatal(FARGS, "Invalid donor %d", i); } - ang *= RAD2DEG; - __ADIST[(int)( ang/abin)]++; - __RDIST[(int)(dist/rbin)]++; - if (!bTwo) + if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET) { - int id, ia; - if ((id = donor_index(&hb->d, grp, i)) == NOTSET) - { - gmx_fatal(FARGS, "Invalid donor %d", i); - } - if ((ia = acceptor_index(&hb->a, ogrp, j)) == NOTSET) - { - gmx_fatal(FARGS, "Invalid acceptor %d", j); - } - resdist = abs(top.atoms.atom[i].resind- - top.atoms.atom[j].resind); - if (resdist >= max_hx) - { - resdist = max_hx-1; - } - __HBDATA->nhx[nframes][resdist]++; + gmx_fatal(FARGS, "Invalid acceptor %d", j); } + resdist = abs(top.atoms.atom[i].resind- + top.atoms.atom[j].resind); + if (resdist >= max_hx) + { + resdist = max_hx-1; + } + __HBDATA->nhx[nframes][resdist]++; } - } - } /* for aj */ - } /* for xjj */ - } /* for yjj */ - } /* for zjj */ - } /* for ai */ - } /* for grp */ - } /* for xi,yi,zi */ - } + + } + } /* for aj */ + } /* for xjj */ + } /* for yjj */ + } /* for zjj */ + } /* for ai */ + } /* for grp */ + } /* for xi,yi,zi */ } - } /* if (bSelected) {...} else */ + } + } /* if (bSelected) {...} else */ - /* Better wait for all threads to finnish using x[] before updating it. */ - k = nframes; + /* Better wait for all threads to finnish using x[] before updating it. */ + k = nframes; #pragma omp barrier #pragma omp critical + { + /* Sum up histograms and counts from p_hb[] into hb */ + if (bOMP) { - /* Sum up histograms and counts from p_hb[] into hb */ - if (bOMP) + hb->nhb[k] += p_hb[threadNr]->nhb[k]; + hb->ndist[k] += p_hb[threadNr]->ndist[k]; + for (j = 0; j < max_hx; j++) { - hb->nhb[k] += p_hb[threadNr]->nhb[k]; - hb->ndist[k] += p_hb[threadNr]->ndist[k]; - for (j = 0; j < max_hx; j++) - { - hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; - } + hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; } } + } - /* Here are a handful of single constructs - * to share the workload a bit. The most - * important one is of course the last one, - * where there's a potential bottleneck in form - * of slow I/O. */ + /* Here are a handful of single constructs + * to share the workload a bit. The most + * important one is of course the last one, + * where there's a potential bottleneck in form + * of slow I/O. */ #pragma omp barrier #pragma omp single - { - analyse_donor_properties(donor_properties, hb, k, t); - } + { + analyse_donor_properties(donor_properties, hb, k, t); + } #pragma omp single + { + if (fpnhb) { - if (fpnhb) - { - do_nhb_dist(fpnhb, hb, t); - } + do_nhb_dist(fpnhb, hb, t); } - } /* if (bNN) {...} else + */ + } #pragma omp single { @@ -4454,11 +3183,8 @@ int gmx_hbond(int argc, char *argv[]) { max_nhb = 0.5*(hb->d.nrd*hb->a.nra); } - /* Added support for -contact below. - * - Erik Marklund, May 29-31, 2006 */ - /* Changed contact code. - * - Erik Marklund, June 29, 2006 */ - if (bHBmap && !bNN) + + if (bHBmap) { if (hb->nrhb == 0) { @@ -4471,8 +3197,6 @@ int gmx_hbond(int argc, char *argv[]) hb->nrhb, bContact ? "contacts" : "hydrogen bonds", hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding"); - /*Control the pHist.*/ - if (bMerge) { merge_hb(hb, bTwo, bContact); @@ -4572,12 +3296,10 @@ int gmx_hbond(int argc, char *argv[]) } xvgrclose(fp); } - if (!bNN) - { - printf("Average number of %s per timeframe %.3f out of %g possible\n", - bContact ? "contacts" : "hbonds", - bContact ? aver_dist : aver_nhb, max_nhb); - } + + printf("Average number of %s per timeframe %.3f out of %g possible\n", + bContact ? "contacts" : "hbonds", + bContact ? aver_dist : aver_nhb, max_nhb); /* Do Autocorrelation etc. */ if (hb->bHBmap) @@ -4596,19 +3318,9 @@ int gmx_hbond(int argc, char *argv[]) { char *gemstring = NULL; - if (bGem || bNN) - { - params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end, - gemBallistic, nBalExp); - if (params == NULL) - { - gmx_fatal(FARGS, "Could not initiate t_gemParams params."); - } - } - gemstring = gmx_strdup(gemType[hb->per->gemtype]); do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump, bMerge, bContact, fit_start, temp, r2cut > 0, oenv, - gemstring, nThreads, NN, bBallistic, bGemFit); + nThreads); } if (opt2bSet("-life", NFILE, fnm)) { @@ -4641,13 +3353,6 @@ int gmx_hbond(int argc, char *argv[]) { if (ISHB(hb->hbmap[id][ia]->history[hh])) { - /* Changed '<' into '<=' in the for-statement below. - * It fixed the previously undiscovered bug that caused - * the last occurance of an hbond/contact to not be - * set in mat.matrix. Have a look at any old -hbm-output - * and you will notice that the last column is allways empty. - * - Erik Marklund May 30, 2006 - */ for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++) { int nn0 = hb->hbmap[id][ia]->n0; @@ -4698,35 +3403,6 @@ int gmx_hbond(int argc, char *argv[]) } } - if (bGem) - { - fprintf(stderr, "There were %i periodic shifts\n", hb->per->nper); - fprintf(stderr, "Freeing pHist for all donors...\n"); - for (i = 0; i < hb->d.nrd; i++) - { - fprintf(stderr, "\r%i", i); - if (hb->per->pHist[i] != NULL) - { - for (j = 0; j < hb->a.nra; j++) - { - clearPshift(&(hb->per->pHist[i][j])); - } - sfree(hb->per->pHist[i]); - } - } - sfree(hb->per->pHist); - sfree(hb->per->p2i); - sfree(hb->per); - fprintf(stderr, "...done.\n"); - } - -#ifdef HAVE_NN_LOOPS - if (bNN) - { - free_hbEmap(hb); - } -#endif - if (hb->bDAnr) { int i, j, nleg; diff --git a/src/gromacs/gmxana/hxprops.c b/src/gromacs/gmxana/hxprops.cpp similarity index 94% rename from src/gromacs/gmxana/hxprops.c rename to src/gromacs/gmxana/hxprops.cpp index 5b1f6c7074..de026aac37 100644 --- a/src/gromacs/gmxana/hxprops.c +++ b/src/gromacs/gmxana/hxprops.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -38,8 +38,10 @@ #include "hxprops.h" -#include -#include +#include +#include + +#include #include "gromacs/legacyheaders/macros.h" #include "gromacs/legacyheaders/typedefs.h" @@ -73,7 +75,6 @@ real ellipticity(int nres, t_bb bb[]) #define NPPW asize(ppw) int i, j; - const real Xi = 1.0; real ell, pp2, phi, psi; ell = 0; @@ -129,7 +130,7 @@ real radius(FILE *fp, int nca, atom_id ca_index[], rvec x[]) fprintf(fp, "\n"); } - return sqrt(dlt/nca); + return std::sqrt(dlt/nca); } static real rot(rvec x1, rvec x2) @@ -137,13 +138,13 @@ static real rot(rvec x1, rvec x2) real phi1, dphi, cp, sp; real xx, yy; - phi1 = atan2(x1[YY], x1[XX]); - cp = cos(phi1); - sp = sin(phi1); + phi1 = std::atan2(x1[YY], x1[XX]); + cp = std::cos(phi1); + sp = std::sin(phi1); xx = cp*x2[XX]+sp*x2[YY]; yy = -sp*x2[XX]+cp*x2[YY]; - dphi = RAD2DEG*atan2(yy, xx); + dphi = RAD2DEG*std::atan2(yy, xx); return dphi; } @@ -341,8 +342,8 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0, r0 = r1 = atom[(*index)[0]].resind; for (i = 1; (i < gnx); i++) { - r0 = min(r0, atom[(*index)[i]].resind); - r1 = max(r1, atom[(*index)[i]].resind); + r0 = std::min(r0, atom[(*index)[i]].resind); + r1 = std::max(r1, atom[(*index)[i]].resind); } rnr = r1-r0+1; fprintf(stderr, "There are %d residues\n", rnr); @@ -356,16 +357,16 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0, { ai = (*index)[i]; ri = atom[ai].resind-r0; - if (strcmp(*(resinfo[ri].name), "PRO") == 0) + if (std::strcmp(*(resinfo[ri].name), "PRO") == 0) { - if (strcmp(*(atomname[ai]), "CD") == 0) + if (std::strcmp(*(atomname[ai]), "CD") == 0) { bb[ri].H = ai; } } for (k = 0; (k < NBB); k++) { - if (strcmp(bb_nm[k], *(atomname[ai])) == 0) + if (std::strcmp(bb_nm[k], *(atomname[ai])) == 0) { j++; break; @@ -427,7 +428,7 @@ t_bb *mkbbind(const char *fn, int *nres, int *nbb, int res0, bb[i].Cprev = bb[i-1].C; bb[i].Nnext = bb[i+1].N; } - rnr = max(0, i1-i0+1); + rnr = std::max(0, i1-i0+1); fprintf(stderr, "There are %d complete backbone residues (from %d to %d)\n", rnr, bb[i0].resno, bb[i1].resno); if (rnr == 0) @@ -462,7 +463,7 @@ real pprms(FILE *fp, int nbb, t_bb bb[]) { if (bb[i].bHelix) { - rms = sqrt(bb[i].pprms2); + rms = std::sqrt(bb[i].pprms2); rmst += rms; rms2 += bb[i].pprms2; fprintf(fp, "%10g ", rms); @@ -470,7 +471,7 @@ real pprms(FILE *fp, int nbb, t_bb bb[]) } } fprintf(fp, "\n"); - rms = sqrt(rms2/n-sqr(rmst/n)); + rms = std::sqrt(rms2/n-sqr(rmst/n)); return rms; } @@ -515,9 +516,9 @@ void calc_hxprops(int nres, t_bb bb[], rvec x[]) bb[i].pprms2 = sqr(bb[i].phi-PHI_AHX)+sqr(bb[i].psi-PSI_AHX); bb[i].jcaha += - 1.4*sin((bb[i].psi+138.0)*DEG2RAD) - - 4.1*cos(2.0*DEG2RAD*(bb[i].psi+138.0)) + - 2.0*cos(2.0*DEG2RAD*(bb[i].phi+30.0)); + 1.4*std::sin((bb[i].psi+138.0)*DEG2RAD) - + 4.1*std::cos(2.0*DEG2RAD*(bb[i].psi+138.0)) + + 2.0*std::cos(2.0*DEG2RAD*(bb[i].phi+30.0)); } } diff --git a/src/gromacs/gmxana/hxprops.h b/src/gromacs/gmxana/hxprops.h index c6242d76b7..9499c5b6bc 100644 --- a/src/gromacs/gmxana/hxprops.h +++ b/src/gromacs/gmxana/hxprops.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -42,6 +42,11 @@ #include "gromacs/legacyheaders/typedefs.h" +#ifdef __cplusplus +extern "C" +{ +#endif + #define PHI_AHX (-55.0) #define PSI_AHX (-45.0) /* Canonical values of the helix phi/psi angles */ @@ -111,4 +116,8 @@ extern void calc_hxprops(int nres, t_bb bb[], rvec x[]); extern void pr_bb(FILE *fp, int nres, t_bb bb[]); +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/gromacs/gmxana/nrama.c b/src/gromacs/gmxana/nrama.cpp similarity index 99% rename from src/gromacs/gmxana/nrama.c rename to src/gromacs/gmxana/nrama.cpp index 165fb75b23..04bfb5ba9a 100644 --- a/src/gromacs/gmxana/nrama.c +++ b/src/gromacs/gmxana/nrama.cpp @@ -38,8 +38,7 @@ #include "nrama.h" -#include -#include +#include #include "gromacs/listed-forces/bonded.h" #include "gromacs/pbcutil/rmpbc.h" @@ -283,7 +282,6 @@ t_topology *init_rama(const output_env_t oenv, const char *infile, const char *topfile, t_xrama *xr, int mult) { t_topology *top; - int ePBC; real t; top = read_top(topfile, &xr->ePBC); diff --git a/src/gromacs/gmxana/nsfactor.c b/src/gromacs/gmxana/nsfactor.cpp similarity index 86% rename from src/gromacs/gmxana/nsfactor.c rename to src/gromacs/gmxana/nsfactor.cpp index 5525de64df..18b682fc4a 100644 --- a/src/gromacs/gmxana/nsfactor.c +++ b/src/gromacs/gmxana/nsfactor.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -38,7 +38,8 @@ #include "config.h" -#include +#include +#include #include "gromacs/fileio/strdb.h" #include "gromacs/math/vec.h" @@ -208,7 +209,7 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( int nthreads; gmx_rng_t *trng = NULL; #endif - gmx_int64_t mc = 0, max; + gmx_int64_t mc = 0, mc_max; gmx_rng_t rng = NULL; /* allocate memory for pr */ @@ -225,8 +226,7 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( rmax = norm(dist); - pr->grn = (int)floor(rmax/pr->binwidth)+1; - rmax = pr->grn*pr->binwidth; + pr->grn = static_cast(std::floor(rmax/pr->binwidth)+1); snew(pr->gr, pr->grn); @@ -235,11 +235,11 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */ if (mcover == -1) { - max = (gmx_int64_t)floor(0.5*0.01*isize*(isize-1)); + mc_max = static_cast(std::floor(0.5*0.01*isize*(isize-1))); } else { - max = (gmx_int64_t)floor(0.5*mcover*isize*(isize-1)); + mc_max = static_cast(std::floor(0.5*mcover*isize*(isize-1))); } rng = gmx_rng_init(seed); #ifdef GMX_OPENMP @@ -256,13 +256,13 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( tid = gmx_omp_get_thread_num(); /* now starting parallel threads */ #pragma omp for - for (mc = 0; mc < max; mc++) + for (mc = 0; mc < mc_max; mc++) { - i = (int)floor(gmx_rng_uniform_real(trng[tid])*isize); - j = (int)floor(gmx_rng_uniform_real(trng[tid])*isize); + i = static_cast(std::floor(gmx_rng_uniform_real(trng[tid])*isize)); + j = static_cast(std::floor(gmx_rng_uniform_real(trng[tid])*isize)); if (i != j) { - tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]]; + tgr[tid][static_cast(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]]; } } } @@ -283,13 +283,13 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( sfree(tgr); sfree(trng); #else - for (mc = 0; mc < max; mc++) + for (mc = 0; mc < mc_max; mc++) { - i = (int)floor(gmx_rng_uniform_real(rng)*isize); - j = (int)floor(gmx_rng_uniform_real(rng)*isize); + i = static_cast(std::floor(gmx_rng_uniform_real(rng)*isize)); + j = static_cast(std::floor(gmx_rng_uniform_real(rng)*isize)); if (i != j) { - pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]]; + pr->gr[static_cast(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]]; } } #endif @@ -314,7 +314,7 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( { for (j = 0; j < i; j++) { - tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]]; + tgr[tid][static_cast(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]]; } } } @@ -337,7 +337,7 @@ gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram ( { for (j = 0; j < i; j++) { - pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]]; + pr->gr[static_cast(std::floor(std::sqrt(distance2(x[index[i]], x[index[j]]))/binwidth))] += gsans->slength[index[i]]*gsans->slength[index[j]]; } } #endif @@ -364,7 +364,7 @@ gmx_static_structurefactor_t *convert_histogram_to_intensity_curve (gmx_radial_d int i, j; /* init data */ snew(sq, 1); - sq->qn = (int)floor((end_q-start_q)/q_step); + sq->qn = static_cast(std::floor((end_q-start_q)/q_step)); snew(sq->q, sq->qn); snew(sq->s, sq->qn); for (i = 0; i < sq->qn; i++) diff --git a/src/gromacs/gmxana/powerspect.c b/src/gromacs/gmxana/powerspect.cpp similarity index 98% rename from src/gromacs/gmxana/powerspect.c rename to src/gromacs/gmxana/powerspect.cpp index 67674959c7..096052508d 100644 --- a/src/gromacs/gmxana/powerspect.c +++ b/src/gromacs/gmxana/powerspect.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2010,2011,2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. diff --git a/src/gromacs/gmxana/pp2shift.c b/src/gromacs/gmxana/pp2shift.cpp similarity index 96% rename from src/gromacs/gmxana/pp2shift.c rename to src/gromacs/gmxana/pp2shift.cpp index fa12f6d0b6..0ca2c3ad7d 100644 --- a/src/gromacs/gmxana/pp2shift.c +++ b/src/gromacs/gmxana/pp2shift.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -36,8 +36,10 @@ */ #include "gmxpre.h" -#include -#include +#include +#include + +#include #include "gromacs/fileio/matio.h" #include "gromacs/gmxana/gstat.h" @@ -131,8 +133,8 @@ static void dump_sd(const char *fn, t_shiftdata *sd) newdata[i][j] = sd->data[i/nfac][j/nfac]; else*/ newdata[i][j] = interpolate(phi, psi, sd); - lo = min(lo, newdata[i][j]); - hi = max(hi, newdata[i][j]); + lo = std::min(lo, newdata[i][j]); + hi = std::max(hi, newdata[i][j]); } } sprintf(buf, "%s.xpm", fn); diff --git a/src/gromacs/gmxana/princ.c b/src/gromacs/gmxana/princ.cpp similarity index 98% rename from src/gromacs/gmxana/princ.c rename to src/gromacs/gmxana/princ.cpp index fb7dcad1d7..c8af2377b3 100644 --- a/src/gromacs/gmxana/princ.c +++ b/src/gromacs/gmxana/princ.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012,2014, by the GROMACS development team, led by + * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -39,6 +39,8 @@ #include "princ.h" +#include + #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/linearalgebra/nrjac.h" @@ -74,7 +76,7 @@ static void ptrans(char *s, real **inten, real d[], real e[]) z = inten[m][3]; n = x*x+y*y+z*z; fprintf(stderr, "%8s %8.3f %8.3f %8.3f, norm:%8.3f, d:%8.3f, e:%8.3f\n", - s, x, y, z, sqrt(n), d[m], e[m]); + s, x, y, z, std::sqrt(n), d[m], e[m]); } fprintf(stderr, "\n"); } diff --git a/src/gromacs/gmxana/sfactor.c b/src/gromacs/gmxana/sfactor.cpp similarity index 93% rename from src/gromacs/gmxana/sfactor.c rename to src/gromacs/gmxana/sfactor.cpp index 3fcb408faf..4542474520 100644 --- a/src/gromacs/gmxana/sfactor.c +++ b/src/gromacs/gmxana/sfactor.cpp @@ -38,6 +38,10 @@ #include "sfactor.h" +#include + +#include + #include "gromacs/fileio/confio.h" #include "gromacs/fileio/strdb.h" #include "gromacs/fileio/trxio.h" @@ -170,9 +174,9 @@ extern void compute_structure_factor (structure_factor_t * sft, matrix box, k_factor[YY] = 2 * M_PI / box[YY][YY]; k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ]; - maxkx = (int) (end_q / k_factor[XX] + 0.5); - maxky = (int) (end_q / k_factor[YY] + 0.5); - maxkz = (int) (end_q / k_factor[ZZ] + 0.5); + maxkx = static_cast(end_q / k_factor[XX] + 0.5); + maxky = static_cast(end_q / k_factor[YY] + 0.5); + maxkz = static_cast(end_q / k_factor[ZZ] + 0.5); snew (counter, sf->n_angles); @@ -185,7 +189,7 @@ extern void compute_structure_factor (structure_factor_t * sft, matrix box, fprintf(stderr, "\n"); for (i = 0; i < maxkx; i++) { - fprintf (stderr, "\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx); + fprintf (stderr, "\rdone %3.1f%% ", (100.0*(i+1))/maxkx); kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) { @@ -195,10 +199,10 @@ extern void compute_structure_factor (structure_factor_t * sft, matrix box, if (i != 0 || j != 0 || k != 0) { kz = k * k_factor[ZZ]; - krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz)); + krr = std::sqrt (sqr (kx) + sqr (ky) + sqr (kz)); if (krr >= start_q && krr <= end_q) { - kr = (int) (krr/sf->ref_k + 0.5); + kr = static_cast(krr/sf->ref_k + 0.5); if (kr < sf->n_angles) { counter[kr]++; /* will be used for the copmutation @@ -210,8 +214,8 @@ extern void compute_structure_factor (structure_factor_t * sft, matrix box, kdotx = kx * redt[p].x[XX] + ky * redt[p].x[YY] + kz * redt[p].x[ZZ]; - tmpSF[i][j][k].re += cos (kdotx) * asf; - tmpSF[i][j][k].im += sin (kdotx) * asf; + tmpSF[i][j][k].re += std::cos(kdotx) * asf; + tmpSF[i][j][k].im += std::sin(kdotx) * asf; } } } @@ -231,10 +235,10 @@ extern void compute_structure_factor (structure_factor_t * sft, matrix box, { ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) { - kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky) - + sqr (kz)); if (krr >= start_q && krr <= end_q) + kz = k * k_factor[ZZ]; krr = std::sqrt (sqr (kx) + sqr (ky) + + sqr (kz)); if (krr >= start_q && krr <= end_q) { - kr = (int) (krr / sf->ref_k + 0.5); + kr = static_cast(krr / sf->ref_k + 0.5); if (kr < sf->n_angles && counter[kr] != 0) { sf->F[group][kr] += @@ -450,13 +454,11 @@ extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX, structure_factor *sf; rvec *xtop; real **sf_table; - int nsftable; matrix box; - double r_tmp; + real r_tmp; gmx_structurefactors_t *gmx_sf; real *a, *b, c; - int success; snew(a, 4); snew(b, 4); @@ -464,7 +466,7 @@ extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX, gmx_sf = gmx_structurefactors_init(fnDAT); - success = gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c); + gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c); snew (sf, 1); sf->energy = energy; @@ -497,12 +499,12 @@ extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX, snew (red, ng); snew (index_atp, ng); - r_tmp = max (box[XX][XX], box[YY][YY]); - r_tmp = (double) max (box[ZZ][ZZ], r_tmp); + r_tmp = std::max(box[XX][XX], box[YY][YY]); + r_tmp = std::max(box[ZZ][ZZ], r_tmp); sf->ref_k = (2.0 * M_PI) / (r_tmp); /* ref_k will be the reference momentum unit */ - sf->n_angles = (int) (end_q / sf->ref_k + 0.5); + sf->n_angles = static_cast(end_q / sf->ref_k + 0.5); snew (sf->F, ng); for (i = 0; i < ng; i++) @@ -575,7 +577,7 @@ extern void save_data (structure_factor_t *sft, const char *file, int ngrps, * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) -> * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2) */ - A = (double) (i * sf->ref_k) / (2.0 * sf->momentum); + A = static_cast(i * sf->ref_k) / (2.0 * sf->momentum); polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A)); sf->F[g][i] *= polarization_factor; } @@ -605,7 +607,7 @@ extern double CMSF (gmx_structurefactors_t *gsf, int type, int nh, double lambda * vectors. See g_sq.h for a short description of CM fit. */ { - int i, success; + int i; double tmp = 0.0, k2; real *a, *b; real c; @@ -632,7 +634,7 @@ extern double CMSF (gmx_structurefactors_t *gsf, int type, int nh, double lambda else { k2 = (sqr (sin_theta) / sqr (10.0 * lambda)); - success = gmx_structurefactors_get_sf(gsf, type, a, b, &c); + gmx_structurefactors_get_sf(gsf, type, a, b, &c); tmp = c; for (i = 0; (i < 4); i++) { @@ -663,7 +665,7 @@ extern real **gmx_structurefactors_table(gmx_structurefactors_t *gsf, real momen snew (sf_table[i], n_angles); for (j = 0; j < n_angles; j++) { - q = ((double) j * ref_k); + q = static_cast(j * ref_k); /* theta is half the angle between incoming and scattered wavevectors. */ sin_theta = q / (2.0 * momentum); @@ -718,7 +720,7 @@ extern real **compute_scattering_factor_table (gmx_structurefactors_t *gsf, stru /* \hbar \omega \lambda = hc = 1239.842 eV * nm */ - sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) / hc); + sf->momentum = (static_cast(2. * 1000.0 * M_PI * sf->energy) / hc); sf->lambda = hc / (1000.0 * sf->energy); fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda); -- 2.22.0