Move generalized born to C++

Change-Id: Ie4498eed9923b7c4df8589cc1abc75d92155abdb

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 dcfee25ed22557193b6a02208b21eedc813b1c26..ebcd9b327d81b9e5256a28c4ad7adb3dff944833 100644 (file)
--- a/src/gromacs/mdlib/genborn.c
+++ b/src/gromacs/mdlib/genborn.cpp
@@ -39,9 +39,12 @@
 
 #include "gromacs/legacyheaders/genborn.h"
 
-#include <math.h>
 #include <string.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #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<int>(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 a18f2e1ac5cbc267a47a4ec97a3bd0b1f1ff1a11..47b14467edf8cce4b91c4fbd27e2802d6c8ae229 100644 (file)
--- 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 <math.h>
+#include <cmath>
+
+#include <algorithm>
 
 #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]);
             }
         }