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) 2011,2014,2015,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 /* This file is completely threadsafe - keep it that way! */
40 #include "hackblock.h"
44 #include "gromacs/gmxpreprocess/notset.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/mdtypes/md_enums.h"
47 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/utility/fatalerror.h"
49 #include "gromacs/utility/smalloc.h"
51 /* these MUST correspond to the enum in hackblock.h */
52 const char *btsNames
[ebtsNR
] = { "bonds", "angles", "dihedrals", "impropers", "exclusions", "cmap" };
53 const int btsNiatoms
[ebtsNR
] = { 2, 3, 4, 4, 2, 5 };
55 static void free_t_bonded(t_rbonded
*rb
)
59 for (i
= 0; i
< MAXATOMLIST
; i
++)
66 static void free_t_bondeds(t_rbondeds
*rbs
)
70 for (i
= 0; i
< rbs
->nb
; i
++)
72 free_t_bonded(&rbs
->b
[i
]);
79 void free_t_restp(int nrtp
, t_restp
**rtp
)
83 for (i
= 0; i
< nrtp
; i
++)
85 sfree((*rtp
)[i
].resname
);
86 sfree((*rtp
)[i
].atom
);
87 for (j
= 0; j
< (*rtp
)[i
].natom
; j
++)
89 sfree(*(*rtp
)[i
].atomname
[j
]);
90 sfree((*rtp
)[i
].atomname
[j
]);
92 sfree((*rtp
)[i
].atomname
);
93 sfree((*rtp
)[i
].cgnr
);
94 for (j
= 0; j
< ebtsNR
; j
++)
96 free_t_bondeds(&(*rtp
)[i
].rb
[j
]);
102 void free_t_hack(int nh
, t_hack
**h
)
106 for (i
= 0; i
< nh
; i
++)
108 sfree((*h
)[i
].oname
);
109 sfree((*h
)[i
].nname
);
111 for (j
= 0; j
< 4; j
++)
120 void free_t_hackblock(int nhb
, t_hackblock
**hb
)
124 for (i
= 0; i
< nhb
; i
++)
126 sfree((*hb
)[i
].name
);
127 free_t_hack((*hb
)[i
].nhack
, &(*hb
)[i
].hack
);
128 for (j
= 0; j
< ebtsNR
; j
++)
130 free_t_bondeds(&(*hb
)[i
].rb
[j
]);
136 void clear_t_hackblock(t_hackblock
*hb
)
144 for (i
= 0; i
< ebtsNR
; i
++)
147 hb
->rb
[i
].b
= nullptr;
151 void clear_t_hack(t_hack
*hack
)
156 hack
->oname
= nullptr;
157 hack
->nname
= nullptr;
158 hack
->atom
= nullptr;
162 for (i
= 0; i
< 4; i
++)
164 hack
->a
[i
] = nullptr;
166 for (i
= 0; i
< DIM
; i
++)
168 hack
->newx
[i
] = NOTSET
;
172 #define safe_strdup(str) ((str != NULL) ? gmx_strdup(str) : NULL)
174 static void copy_t_rbonded(t_rbonded
*s
, t_rbonded
*d
)
178 for (i
= 0; i
< MAXATOMLIST
; i
++)
180 d
->a
[i
] = safe_strdup(s
->a
[i
]);
182 d
->s
= safe_strdup(s
->s
);
186 static gmx_bool
contains_char(t_rbonded
*s
, char c
)
192 for (i
= 0; i
< MAXATOMLIST
; i
++)
194 if (s
->a
[i
] && s
->a
[i
][0] == c
)
204 rbonded_find_atoms_in_list(t_rbonded
*b
, t_rbonded blist
[], int nlist
, int natoms
)
210 for (i
= 0; i
< nlist
&& foundPos
< 0; i
++)
213 for (k
= 0; k
< natoms
&& atomsMatch
; k
++)
215 atomsMatch
= atomsMatch
&& !strcmp(b
->a
[k
], blist
[i
].a
[k
]);
217 /* Try reverse if forward match did not work */
221 for (k
= 0; k
< natoms
&& atomsMatch
; k
++)
223 atomsMatch
= atomsMatch
&& !strcmp(b
->a
[k
], blist
[i
].a
[natoms
-1-k
]);
229 /* If all the atoms AND all the parameters match, it is likely that
230 * the user made a copy-and-paste mistake (since it would be much cheaper
231 * to just bump the force constant 2x if you really want it twice).
232 * Since we only have the unparsed string here we can only detect
233 * EXACT matches (including identical whitespace).
235 if (!strcmp(b
->s
, blist
[i
].s
))
237 gmx_warning("Duplicate line found in or between hackblock and rtp entries");
244 gmx_bool
merge_t_bondeds(t_rbondeds s
[], t_rbondeds d
[], gmx_bool bMin
, gmx_bool bPlus
)
247 gmx_bool bBondsRemoved
;
248 int nbHackblockStart
;
251 bBondsRemoved
= FALSE
;
252 for (i
= 0; i
< ebtsNR
; i
++)
256 /* Record how many bonds we have in the destination when we start.
258 * If an entry is present in the hackblock (destination), we will
259 * not add the one from the main rtp, since the point is for hackblocks
260 * to overwrite it. However, if there is no hackblock entry we do
261 * allow multiple main rtp entries since some forcefield insist on that.
263 * We accomplish this by checking the position we find an entry in,
264 * rather than merely checking whether it exists at all.
265 * If that index is larger than the original (hackblock) destination
266 * size, it was added from the main rtp, and then we will allow more
267 * such entries. In contrast, if the entry found has a lower index
268 * it is a hackblock entry meant to override the main rtp, and then
269 * we don't add the main rtp one.
271 nbHackblockStart
= d
[i
].nb
;
274 srenew(d
[i
].b
, d
[i
].nb
+ s
[i
].nb
);
275 for (j
= 0; j
< s
[i
].nb
; j
++)
277 /* Check if this bonded string already exists before adding.
278 * We are merging from the main RTP to the hackblocks, so this
279 * will mean the hackblocks overwrite the man RTP, as intended.
281 index
= rbonded_find_atoms_in_list(&s
[i
].b
[j
], d
[i
].b
, d
[i
].nb
, btsNiatoms
[i
]);
282 /* - If we did not find this interaction at all, the index will be -1,
283 * and then we should definitely add it to the merged hackblock and rtp.
285 * Alternatively, if it was found, index will be >=0.
286 * - In case this index is lower than the original number of entries,
287 * it is already present as a *hackblock* entry, and those should
288 * always override whatever we have listed in the RTP. Thus, we
289 * should just keep that one and not add anything from the RTP.
290 * - Finally, if it was found, but with an index higher than
291 * the original number of entries, it comes from the RTP rather
292 * than hackblock, and then we must have added it ourselves
293 * in a previous iteration. In that case it is a matter of
294 * several entries for the same sequence of atoms, and we allow
295 * that in the RTP. In this case we should simply copy all of
296 * them, including this one.
298 if (index
< 0 || index
>= nbHackblockStart
)
300 if (!(bMin
&& contains_char(&s
[i
].b
[j
], '-'))
301 && !(bPlus
&& contains_char(&s
[i
].b
[j
], '+')))
303 copy_t_rbonded(&s
[i
].b
[j
], &d
[i
].b
[ d
[i
].nb
]);
306 else if (i
== ebtsBONDS
)
308 bBondsRemoved
= TRUE
;
313 /* This is the common case where a hackblock entry simply
314 * overrides the RTP, so we cannot warn here.
320 return bBondsRemoved
;
323 void copy_t_restp(t_restp
*s
, t_restp
*d
)
328 d
->resname
= safe_strdup(s
->resname
);
329 snew(d
->atom
, s
->natom
);
330 for (i
= 0; i
< s
->natom
; i
++)
332 d
->atom
[i
] = s
->atom
[i
];
334 snew(d
->atomname
, s
->natom
);
335 for (i
= 0; i
< s
->natom
; i
++)
337 snew(d
->atomname
[i
], 1);
338 *d
->atomname
[i
] = safe_strdup(*s
->atomname
[i
]);
340 snew(d
->cgnr
, s
->natom
);
341 for (i
= 0; i
< s
->natom
; i
++)
343 d
->cgnr
[i
] = s
->cgnr
[i
];
345 for (i
= 0; i
< ebtsNR
; i
++)
347 d
->rb
[i
].type
= s
->rb
[i
].type
;
349 d
->rb
[i
].b
= nullptr;
351 merge_t_bondeds(s
->rb
, d
->rb
, FALSE
, FALSE
);
354 void copy_t_hack(t_hack
*s
, t_hack
*d
)
359 d
->oname
= safe_strdup(s
->oname
);
360 d
->nname
= safe_strdup(s
->nname
);
364 *(d
->atom
) = *(s
->atom
);
370 for (i
= 0; i
< 4; i
++)
372 d
->a
[i
] = safe_strdup(s
->a
[i
]);
374 copy_rvec(s
->newx
, d
->newx
);
377 void merge_hacks_lo(int ns
, t_hack
*s
, int *nd
, t_hack
**d
)
383 srenew(*d
, *nd
+ ns
);
384 for (i
= 0; i
< ns
; i
++)
386 copy_t_hack(&s
[i
], &(*d
)[*nd
+ i
]);
392 void merge_hacks(t_hackblock
*s
, t_hackblock
*d
)
394 merge_hacks_lo(s
->nhack
, s
->hack
, &d
->nhack
, &d
->hack
);
397 void merge_t_hackblock(t_hackblock
*s
, t_hackblock
*d
)
400 merge_t_bondeds(s
->rb
, d
->rb
, FALSE
, FALSE
);
403 void copy_t_hackblock(t_hackblock
*s
, t_hackblock
*d
)
408 d
->name
= safe_strdup(s
->name
);
411 for (i
= 0; i
< ebtsNR
; i
++)
414 d
->rb
[i
].b
= nullptr;
416 merge_t_hackblock(s
, d
);
421 void dump_hb(FILE *out
, int nres
, t_hackblock hb
[])
425 #define SS(s) (s) ? (s) : "-"
426 #define SA(s) (s) ? "+" : ""
427 fprintf(out
, "t_hackblock\n");
428 for (i
= 0; i
< nres
; i
++)
430 fprintf(out
, "%3d %4s %2d %2d\n",
431 i
, SS(hb
[i
].name
), hb
[i
].nhack
, hb
[i
].maxhack
);
434 for (j
= 0; j
< hb
[i
].nhack
; j
++)
436 fprintf(out
, "%d: %d %4s %4s %1s %2d %d %4s %4s %4s %4s\n",
438 SS(hb
[i
].hack
[j
].oname
), SS(hb
[i
].hack
[j
].nname
),
439 SA(hb
[i
].hack
[j
].atom
), hb
[i
].hack
[j
].tp
, hb
[i
].hack
[j
].cgnr
,
440 SS(hb
[i
].hack
[j
].ai()), SS(hb
[i
].hack
[j
].aj()),
441 SS(hb
[i
].hack
[j
].ak()), SS(hb
[i
].hack
[j
].al()) );
444 for (j
= 0; j
< ebtsNR
; j
++)
448 fprintf(out
, " %c %d:", btsNames
[j
][0], hb
[i
].rb
[j
].nb
);
449 for (k
= 0; k
< hb
[i
].rb
[j
].nb
; k
++)
452 for (l
= 0; l
< btsNiatoms
[j
]; l
++)
454 fprintf(out
, " %s", hb
[i
].rb
[j
].b
[k
].a
[l
]);
456 fprintf(out
, " %s]", SS(hb
[i
].rb
[j
].b
[k
].s
));
458 fprintf(out
, " Entry matched: %s\n", yesno_names
[hb
[i
].rb
[j
].b
[k
].match
]);