2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /* This file is completely threadsafe - keep it that way! */
48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
49 #include "gromacs/gmxpreprocess/grompp_impl.h"
50 #include "gromacs/gmxpreprocess/toputil.h"
51 #include "gromacs/topology/ifunc.h"
52 #include "gromacs/utility/cstringutil.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/smalloc.h"
63 static t_2morse
* read_dissociation_energies(int* n2morse
)
67 const char* fn
= "edissoc.dat";
68 t_2morse
* t2m
= nullptr;
69 int maxn2m
= 0, n2m
= 0;
72 /* Open the file with dissociation energies */
73 gmx::FilePtr fp
= gmx::openLibraryFile(fn
);
76 /* Try and read two atom names and an energy term from it */
77 nread
= fscanf(fp
.get(), "%s%s%lf", ai
, aj
, &e_diss
);
80 /* If we got three terms, it probably was OK, no further checking */
83 /* Increase memory for 16 at once, some mallocs are stupid */
88 t2m
[n2m
].ai
= gmx_strdup(ai
);
89 t2m
[n2m
].aj
= gmx_strdup(aj
);
90 t2m
[n2m
].e_diss
= e_diss
;
91 /* Increment counter */
94 /* If we did not read three items, quit reading */
97 /* Set the return values */
103 static int nequal(const char* a1
, const char* a2
)
107 /* Count the number of (case insensitive) characters that are equal in
108 * two strings. If they are equally long their respective null characters are
111 for (i
= 0; (a1
[i
] != '\0') && (a2
[i
] != '\0'); i
++)
113 if (toupper(a1
[i
]) != toupper(a2
[i
]))
118 if ((a1
[i
] == '\0') && (a2
[i
] == '\0'))
126 static real
search_e_diss(int n2m
, t_2morse t2m
[], const char* ai
, const char* aj
)
130 int nii
, njj
, nbstii
= 0, nbstjj
= 0;
133 /* Do a best match search for dissociation energies */
134 for (i
= 0; (i
< n2m
); i
++)
136 /* Check for a perfect match */
137 if (((gmx_strcasecmp(t2m
[i
].ai
, ai
) == 0) && (gmx_strcasecmp(t2m
[i
].aj
, aj
) == 0))
138 || ((gmx_strcasecmp(t2m
[i
].aj
, ai
) == 0) && (gmx_strcasecmp(t2m
[i
].ai
, aj
) == 0)))
145 /* Otherwise count the number of equal characters in the strings ai and aj
146 * and the ones from the file
148 nii
= nequal(t2m
[i
].ai
, ai
);
149 njj
= nequal(t2m
[i
].aj
, aj
);
150 if (((nii
> nbstii
) && (njj
>= nbstjj
)) || ((nii
>= nbstii
) && (njj
> nbstjj
)))
152 if ((nii
> 0) && (njj
> 0))
161 /* Swap ai and aj (at least in counting the number of equal chars) */
162 nii
= nequal(t2m
[i
].ai
, aj
);
163 njj
= nequal(t2m
[i
].aj
, ai
);
164 if (((nii
> nbstii
) && (njj
>= nbstjj
)) || ((nii
>= nbstii
) && (njj
> nbstjj
)))
166 if ((nii
> 0) && (njj
> 0))
176 /* Return the dissocation energy corresponding to the best match, if we have
185 return t2m
[ibest
].e_diss
;
189 void convert_harmonics(gmx::ArrayRef
<MoleculeInformation
> mols
, PreprocessingAtomTypes
* atype
)
194 /* First get the data */
195 t2m
= read_dissociation_energies(&n2m
);
198 fprintf(stderr
, "No dissocation energies read\n");
202 /* For all the molecule types */
204 for (auto& mol
: mols
)
206 /* Check how many morse and harmonic BONDSs there are, increase size of
207 * morse with the number of harmonics
209 for (int bb
= 0; (bb
< F_NRE
); bb
++)
211 if ((interaction_function
[bb
].flags
& IF_BTYPE
) && (bb
!= F_MORSE
))
213 int nrharm
= mol
.interactions
[bb
].size();
215 /* Now loop over the harmonics, trying to convert them */
216 for (auto harmonic
= mol
.interactions
[bb
].interactionTypes
.begin();
217 harmonic
!= mol
.interactions
[bb
].interactionTypes
.end();)
219 int ni
= harmonic
->ai();
220 int nj
= harmonic
->aj();
221 real edis
= search_e_diss(n2m
, t2m
,
222 atype
->atomNameFromAtomType(mol
.atoms
.atom
[ni
].type
),
223 atype
->atomNameFromAtomType(mol
.atoms
.atom
[nj
].type
));
226 real b0
= harmonic
->c0();
227 real kb
= harmonic
->c1();
228 real beta
= std::sqrt(kb
/ (2 * edis
));
229 std::vector
<int> atoms
= { ni
, nj
};
230 std::vector
<real
> forceParam
= { b0
, edis
, beta
};
231 mol
.interactions
[F_MORSE
].interactionTypes
.emplace_back(
232 InteractionOfType(atoms
, forceParam
));
233 harmonic
= mol
.interactions
[bb
].interactionTypes
.erase(harmonic
);
241 int newHarmonics
= mol
.interactions
[bb
].size();
242 fprintf(stderr
, "Converted %d out of %d %s to morse bonds for mol %d\n",
243 nrharm
- newHarmonics
, nrharm
, interaction_function
[bb
].name
, i
);