From a751ebf72b76e72acc733b63faaa0215d660e7de Mon Sep 17 00:00:00 2001 From: Erik Lindahl Date: Wed, 8 Jul 2015 12:57:35 +0200 Subject: [PATCH] Move generalized born to C++ Change-Id: Ie4498eed9923b7c4df8589cc1abc75d92155abdb --- src/gromacs/mdlib/{genborn.c => genborn.cpp} | 94 +++++++------------ ...enborn_allvsall.c => genborn_allvsall.cpp} | 32 +++---- 2 files changed, 46 insertions(+), 80 deletions(-) rename src/gromacs/mdlib/{genborn.c => genborn.cpp} (95%) rename src/gromacs/mdlib/{genborn_allvsall.c => genborn_allvsall.cpp} (97%) diff --git a/src/gromacs/mdlib/genborn.c b/src/gromacs/mdlib/genborn.cpp similarity index 95% rename from src/gromacs/mdlib/genborn.c rename to src/gromacs/mdlib/genborn.cpp index dcfee25ed2..ebcd9b327d 100644 --- a/src/gromacs/mdlib/genborn.c +++ b/src/gromacs/mdlib/genborn.cpp @@ -39,9 +39,12 @@ #include "gromacs/legacyheaders/genborn.h" -#include #include +#include + +#include + #include "gromacs/domdec/domdec.h" #include "gromacs/fileio/pdbio.h" #include "gromacs/legacyheaders/names.h" @@ -119,12 +122,8 @@ static int init_gb_still(const t_atomtypes *atype, t_idef *idef, t_atoms *atoms, gmx_genborn_t *born, int natoms) { - int i, j, i1, i2, k, m, nbond, nang, ia, ib, ic, id, nb, idx, idx2, at; - int iam, ibm; - int at0, at1; - real length, angle; - real r, ri, rj, ri2, ri3, rj2, r2, r3, r4, rk, ratio, term, h, doffset; - real p1, p2, p3, factor, cosine, rab, rbc; + int i, j, m, ia, ib; + real r, ri, rj, ri2, rj2, r3, r4, ratio, term, h, doffset; real *vsol; real *gp; @@ -133,9 +132,6 @@ static int init_gb_still(const t_atomtypes *atype, t_idef *idef, t_atoms *atoms, snew(gp, natoms); snew(born->gpol_still_work, natoms+3); - at0 = 0; - at1 = natoms; - doffset = born->gb_doffset; for (i = 0; i < natoms; i++) @@ -165,7 +161,6 @@ static int init_gb_still(const t_atomtypes *atype, t_idef *idef, t_atoms *atoms, rj = atype->gb_radius[atoms->atom[ib].type]; ri2 = ri*ri; - ri3 = ri2*ri; rj2 = rj*rj; ratio = (rj2-ri2-r*r)/(2*ri*r); @@ -237,8 +232,8 @@ int init_gb(gmx_genborn_t **p_born, t_forcerec *fr, const t_inputrec *ir, const gmx_mtop_t *mtop, int gb_algorithm) { - int i, j, m, ai, aj, jj, natoms, nalloc; - real rai, sk, p, doffset; + int i, jj, natoms; + real rai, sk, doffset; t_atoms atoms; gmx_genborn_t *born; @@ -363,10 +358,10 @@ calc_gb_rad_still(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md) { - int i, k, n, nj0, nj1, ai, aj, type; + int i, k, n, nj0, nj1, ai, aj; int shift; real shX, shY, shZ; - real gpi, dr, dr2, dr4, idr4, rvdw, ratio, ccf, theta, term, rai, raj; + real gpi, dr2, idr4, rvdw, ratio, ccf, theta, term, rai, raj; real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11; real rinv, idr2, idr6, vaj, dccf, cosq, sinq, prod, gpi2; real factor; @@ -491,14 +486,14 @@ calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md) { - int i, k, n, ai, aj, nj0, nj1, at0, at1; + int i, k, n, ai, aj, nj0, nj1; int shift; real shX, shY, shZ; - real rai, raj, gpi, dr2, dr, sk, sk_ai, sk2, sk2_ai, lij, uij, diff2, tmp, sum_ai; + real rai, raj, dr2, dr, sk, sk_ai, sk2, sk2_ai, lij, uij, diff2, tmp, sum_ai; real rad, min_rad, rinv, rai_inv; real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11; real lij2, uij2, lij3, uij3, t1, t2, t3; - real lij_inv, dlij, duij, sk2_rinv, prod, log_term; + real lij_inv, dlij, sk2_rinv, prod, log_term; real doffset, raj_inv, dadx_val; real *gb_radius; @@ -512,7 +507,6 @@ calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, /* Keep the compiler happy */ n = 0; - prod = 0; for (i = 0; i < nl->nri; i++) { @@ -585,7 +579,7 @@ calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, sk2_rinv = sk2*rinv; prod = 0.25*sk2_rinv; - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -641,7 +635,7 @@ calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, /* log_term = table_log(uij*lij_inv,born->log_table, LOG_TABLE_ACCURACY); */ - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -686,7 +680,7 @@ calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, min_rad = rai + doffset; rad = 1.0/sum_ai; - born->bRad[i] = rad > min_rad ? rad : min_rad; + born->bRad[i] = std::max(rad, min_rad); fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]); } } @@ -706,21 +700,19 @@ static int calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md) { - int i, k, ai, aj, nj0, nj1, n, at0, at1; + int i, k, ai, aj, nj0, nj1, n; int shift; real shX, shY, shZ; - real rai, raj, gpi, dr2, dr, sk, sk2, lij, uij, diff2, tmp, sum_ai; - real rad, min_rad, sum_ai2, sum_ai3, tsum, tchain, rinv, rai_inv, lij_inv, rai_inv2; + real rai, raj, dr2, dr, sk, sk2, lij, uij, diff2, tmp, sum_ai; + real sum_ai2, sum_ai3, tsum, tchain, rinv, rai_inv, lij_inv, rai_inv2; real log_term, prod, sk2_rinv, sk_ai, sk2_ai; real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11; - real lij2, uij2, lij3, uij3, dlij, duij, t1, t2, t3; + real lij2, uij2, lij3, uij3, dlij, t1, t2, t3; real doffset, raj_inv, dadx_val; real *gb_radius; /* Keep the compiler happy */ n = 0; - prod = 0; - raj = 0; doffset = born->gb_doffset; gb_radius = born->gb_radius; @@ -801,7 +793,7 @@ calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, sk2_rinv = sk2*rinv; prod = 0.25*sk2_rinv; - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -853,7 +845,7 @@ calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, prod = 0.25 * sk2_rinv; /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */ - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -928,7 +920,6 @@ calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top, int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir, gmx_localtop_t *top, rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md, t_nrnb *nrnb) { - real *p; int cnt; int ndadx; @@ -1050,14 +1041,14 @@ real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtab real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r, real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph) { - int i, j, n0, m, nnn, type, ai, aj; + int i, j, n0, m, nnn, ai, aj; int ki; real isai, isaj; real r, rsq11; real rinv11, iq; real isaprod, qq, gbscale, gbtabscale, Y, F, Geps, Heps2, Fp, VV, FF, rt, eps, eps2; - real vgb, fgb, vcoul, fijC, dvdatmp, fscal, dvdaj; + real vgb, fgb, fijC, dvdatmp, fscal; real vctot; rvec dx; @@ -1097,7 +1088,7 @@ real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtab gbscale = isaprod*gbtabscale; r = rsq11*rinv11; rt = r*gbscale; - n0 = rt; + n0 = static_cast(rt); eps = rt-n0; eps2 = eps*eps; nnn = 4*n0; @@ -1186,12 +1177,9 @@ real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr, int natoms, gmx_genborn_t * real *dvda, t_mdatoms *md) { int ai, i, at0, at1; - real e, es, rai, rbi, term, probe, tmp, factor; + real e, es, rai, term, probe, tmp, factor; real rbi_inv, rbi_inv2; - /* To keep the compiler happy */ - factor = 0; - if (DOMAINDECOMP(cr)) { at0 = 0; @@ -1205,19 +1193,6 @@ real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr, int natoms, gmx_genborn_t * /* factor is the surface tension */ factor = born->sa_surface_tension; - /* - - // The surface tension factor is 0.0049 for Still model, 0.0054 for HCT/OBC - if(gb_algorithm==egbSTILL) - { - factor=0.0049*100*CAL2JOULE; - } - else - { - factor=0.0054*100*CAL2JOULE; - } - */ - /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */ es = 0; probe = 0.14; @@ -1251,11 +1226,10 @@ real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[ int i, k, n, ai, aj, nj0, nj1, n0, n1; int shift; real shX, shY, shZ; - real fgb, fij, rb2, rbi, fix1, fiy1, fiz1; - real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11, rsq11; - real rinv11, tx, ty, tz, rbai, rbaj, fgb_ai; + real fgb, rbi, fix1, fiy1, fiz1; + real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11; + real tx, ty, tz, rbai, rbaj, fgb_ai; real *rb; - volatile int idx; n = 0; rb = born->work; @@ -1366,9 +1340,7 @@ calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb, const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd) { - real v = 0; int cnt; - int i; /* PBC or not? */ const t_pbc *pbc_null; @@ -1468,18 +1440,17 @@ static void add_bondeds_to_gblist(t_ilist *il, gmx_bool bMolPBC, t_pbc *pbc, t_graph *g, rvec *x, struct gbtmpnbls *nls) { - int ind, j, ai, aj, shift, found; + int ind, j, ai, aj, found; rvec dx; ivec dt; gbtmpnbl_t *list; - shift = CENTRAL; for (ind = 0; ind < il->nr; ind += 3) { ai = il->iatoms[ind+1]; aj = il->iatoms[ind+2]; - shift = CENTRAL; + int shift = CENTRAL; if (g != NULL) { rvec_sub(x[ai], x[aj], dx); @@ -1531,8 +1502,7 @@ int make_gb_nblist(t_commrec *cr, int gb_algorithm, rvec x[], matrix box, t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born) { - int i, l, ii, j, k, n, nj0, nj1, ai, aj, at0, at1, found, shift, s; - int apa; + int i, j, k, n, nj0, nj1, ai, shift, s; t_nblist *nblist; t_pbc pbc; diff --git a/src/gromacs/mdlib/genborn_allvsall.c b/src/gromacs/mdlib/genborn_allvsall.cpp similarity index 97% rename from src/gromacs/mdlib/genborn_allvsall.c rename to src/gromacs/mdlib/genborn_allvsall.cpp index a18f2e1ac5..47b14467ed 100644 --- a/src/gromacs/mdlib/genborn_allvsall.c +++ b/src/gromacs/mdlib/genborn_allvsall.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2009, The GROMACS Development Team. - * Copyright (c) 2010,2014, by the GROMACS development team, led by + * Copyright (c) 2010,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,9 @@ #include "genborn_allvsall.h" -#include +#include + +#include #include "gromacs/legacyheaders/genborn.h" #include "gromacs/legacyheaders/network.h" @@ -115,12 +117,10 @@ setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t * aadata, gmx_bool bInclude13, gmx_bool bInclude14) { - int i, j, k, tp; + int i, j, k; int a1, a2; - int nj0, nj1; int max_offset; int max_excl_offset; - int nj; /* This routine can appear to be a bit complex, but it is mostly book-keeping. * To enable the fast all-vs-all kernel we need to be able to stream through all coordinates @@ -168,7 +168,7 @@ setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t * aadata, } if (k > 0 && k <= max_offset) { - max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset; + max_excl_offset = std::max(k, max_excl_offset); } } } @@ -194,7 +194,7 @@ setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t * aadata, } if (k > 0 && k <= max_offset) { - max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset; + max_excl_offset = std::max(k, max_excl_offset); } } } @@ -220,11 +220,11 @@ setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t * aadata, } if (k > 0 && k <= max_offset) { - max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset; + max_excl_offset = std::max(k, max_excl_offset); } } } - max_excl_offset = (max_offset < max_excl_offset) ? max_offset : max_excl_offset; + max_excl_offset = std::min(max_offset, max_excl_offset); aadata->jindex_gb[3*i+1] = i+1+max_excl_offset; @@ -329,9 +329,7 @@ genborn_allvsall_setup(gmx_allvsallgb2_data_t ** p_aadata, gmx_bool bInclude13, gmx_bool bInclude14) { - int i, j, idx; gmx_allvsallgb2_data_t *aadata; - real *p; snew(aadata, 1); *p_aadata = aadata; @@ -584,8 +582,6 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, ni1 = mdatoms->homenr; n = 0; - prod = 0; - raj = 0; doffset = born->gb_doffset; aadata = *((gmx_allvsallgb2_data_t **)work); @@ -679,7 +675,7 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, sk2_rinv = sk2*rinv; prod = 0.25*sk2_rinv; - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */ tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -730,7 +726,7 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, prod = 0.25 * sk2_rinv; /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */ - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -806,7 +802,7 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, sk2_rinv = sk2*rinv; prod = 0.25*sk2_rinv; - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */ tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -857,7 +853,7 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, prod = 0.25 * sk2_rinv; /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */ - log_term = log(uij*lij_inv); + log_term = std::log(uij*lij_inv); tmp = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2); @@ -898,7 +894,7 @@ genborn_allvsall_calc_hct_obc_radii(t_forcerec * fr, min_rad = rai + born->gb_doffset; rad = 1.0/sum_ai; - born->bRad[i] = rad > min_rad ? rad : min_rad; + born->bRad[i] = std::max(rad, min_rad); fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]); } } -- 2.22.0