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, 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.
49 #include <sys/types.h>
51 #include "gromacs/fileio/gmxfio.h"
52 #include "gromacs/fileio/txtdump.h"
53 #include "gromacs/gmxlib/warninp.h"
54 #include "gromacs/gmxpreprocess/gmxcpp.h"
55 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
56 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
57 #include "gromacs/gmxpreprocess/grompp-impl.h"
58 #include "gromacs/gmxpreprocess/topdirs.h"
59 #include "gromacs/gmxpreprocess/toppush.h"
60 #include "gromacs/gmxpreprocess/topshake.h"
61 #include "gromacs/gmxpreprocess/toputil.h"
62 #include "gromacs/gmxpreprocess/vsite_parm.h"
63 #include "gromacs/legacyheaders/names.h"
64 #include "gromacs/legacyheaders/types/ifunc.h"
65 #include "gromacs/math/units.h"
66 #include "gromacs/math/utilities.h"
67 #include "gromacs/mdlib/genborn.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/symtab.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/futil.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/smalloc.h"
78 #define OPENDIR '[' /* starting sign for directive */
79 #define CLOSEDIR ']' /* ending sign for directive */
81 static void free_nbparam(t_nbparam
**param
, int nr
)
85 for (i
= 0; i
< nr
; i
++)
92 static int copy_nbparams(t_nbparam
**param
, int ftype
, t_params
*plist
, int nr
)
100 for (i
= 0; i
< nr
; i
++)
102 for (j
= 0; j
<= i
; j
++)
104 if (param
[i
][j
].bSet
)
106 for (f
= 0; f
< nrfp
; f
++)
108 plist
->param
[nr
*i
+j
].c
[f
] = param
[i
][j
].c
[f
];
109 plist
->param
[nr
*j
+i
].c
[f
] = param
[i
][j
].c
[f
];
119 static void gen_pairs(t_params
*nbs
, t_params
*pairs
, real fudge
, int comb
)
121 int i
, j
, ntp
, nrfp
, nrfpA
, nrfpB
, nnn
;
124 nnn
= static_cast<int>(std::sqrt(static_cast<double>(ntp
)));
125 GMX_ASSERT(nnn
* nnn
== ntp
, "Number of pairs of generated non-bonded parameters should be a perfect square");
127 nrfpA
= interaction_function
[F_LJ14
].nrfpA
;
128 nrfpB
= interaction_function
[F_LJ14
].nrfpB
;
131 if ((nrfp
!= nrfpA
) || (nrfpA
!= nrfpB
))
133 gmx_incons("Number of force parameters in gen_pairs wrong");
136 fprintf(stderr
, "Generating 1-4 interactions: fudge = %g\n", fudge
);
139 fprintf(debug
, "Fudge factor for 1-4 interactions: %g\n", fudge
);
140 fprintf(debug
, "Holy Cow! there are %d types\n", ntp
);
142 snew(pairs
->param
, pairs
->nr
);
143 for (i
= 0; (i
< ntp
); i
++)
146 pairs
->param
[i
].a
[0] = i
/ nnn
;
147 pairs
->param
[i
].a
[1] = i
% nnn
;
148 /* Copy normal and FEP parameters and multiply by fudge factor */
152 for (j
= 0; (j
< nrfp
); j
++)
154 /* If we are using sigma/epsilon values, only the epsilon values
155 * should be scaled, but not sigma.
156 * The sigma values have even indices 0,2, etc.
158 if ((comb
== eCOMB_ARITHMETIC
|| comb
== eCOMB_GEOM_SIG_EPS
) && (j
%2 == 0))
167 pairs
->param
[i
].c
[j
] = scaling
*nbs
->param
[i
].c
[j
];
168 pairs
->param
[i
].c
[nrfp
+j
] = scaling
*nbs
->param
[i
].c
[j
];
173 double check_mol(gmx_mtop_t
*mtop
, warninp_t wi
)
176 int i
, mb
, nmol
, ri
, pt
;
181 /* Check mass and charge */
184 for (mb
= 0; mb
< mtop
->nmoltype
; mb
++)
186 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
187 nmol
= mtop
->molblock
[mb
].nmol
;
188 for (i
= 0; (i
< atoms
->nr
); i
++)
190 q
+= nmol
*atoms
->atom
[i
].q
;
191 m
= atoms
->atom
[i
].m
;
192 mB
= atoms
->atom
[i
].mB
;
193 pt
= atoms
->atom
[i
].ptype
;
194 /* If the particle is an atom or a nucleus it must have a mass,
195 * else, if it is a shell, a vsite or a bondshell it can have mass zero
197 if (((m
<= 0.0) || (mB
<= 0.0)) && ((pt
== eptAtom
) || (pt
== eptNucleus
)))
199 ri
= atoms
->atom
[i
].resind
;
200 sprintf(buf
, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
201 *(atoms
->atomname
[i
]),
202 *(atoms
->resinfo
[ri
].name
),
203 atoms
->resinfo
[ri
].nr
,
205 warning_error(wi
, buf
);
208 if (((m
!= 0) || (mB
!= 0)) && (pt
== eptVSite
))
210 ri
= atoms
->atom
[i
].resind
;
211 sprintf(buf
, "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state B)\n"
212 " Check your topology.\n",
213 *(atoms
->atomname
[i
]),
214 *(atoms
->resinfo
[ri
].name
),
215 atoms
->resinfo
[ri
].nr
,
217 warning_error(wi
, buf
);
218 /* The following statements make LINCS break! */
219 /* atoms->atom[i].m=0; */
226 static void sum_q(t_atoms
*atoms
, int n
, double *qt
, double *qBt
)
234 for (i
= 0; i
< atoms
->nr
; i
++)
236 qmolA
+= atoms
->atom
[i
].q
;
237 qmolB
+= atoms
->atom
[i
].qB
;
239 /* Unfortunately an absolute comparison,
240 * but this avoids unnecessary warnings and gmx-users mails.
242 if (fabs(qmolA
) >= 1e-6 || fabs(qmolB
) >= 1e-6)
249 static void get_nbparm(char *nb_str
, char *comb_str
, int *nb
, int *comb
,
253 char warn_buf
[STRLEN
];
256 for (i
= 1; (i
< eNBF_NR
); i
++)
258 if (gmx_strcasecmp(nb_str
, enbf_names
[i
]) == 0)
265 *nb
= strtol(nb_str
, NULL
, 10);
267 if ((*nb
< 1) || (*nb
>= eNBF_NR
))
269 sprintf(warn_buf
, "Invalid nonbond function selector '%s' using %s",
270 nb_str
, enbf_names
[1]);
271 warning_error(wi
, warn_buf
);
275 for (i
= 1; (i
< eCOMB_NR
); i
++)
277 if (gmx_strcasecmp(comb_str
, ecomb_names
[i
]) == 0)
284 *comb
= strtol(comb_str
, NULL
, 10);
286 if ((*comb
< 1) || (*comb
>= eCOMB_NR
))
288 sprintf(warn_buf
, "Invalid combination rule selector '%s' using %s",
289 comb_str
, ecomb_names
[1]);
290 warning_error(wi
, warn_buf
);
295 static char ** cpp_opts(const char *define
, const char *include
,
300 const char *cppadds
[2];
301 char **cppopts
= NULL
;
302 const char *option
[2] = { "-D", "-I" };
303 const char *nopt
[2] = { "define", "include" };
307 char warn_buf
[STRLEN
];
310 cppadds
[1] = include
;
311 for (n
= 0; (n
< 2); n
++)
318 while ((*ptr
!= '\0') && isspace(*ptr
))
323 while ((*rptr
!= '\0') && !isspace(*rptr
))
331 strncpy(buf
, ptr
, len
);
332 if (strstr(ptr
, option
[n
]) != ptr
)
334 set_warning_line(wi
, "mdp file", -1);
335 sprintf(warn_buf
, "Malformed %s option %s", nopt
[n
], buf
);
336 warning(wi
, warn_buf
);
340 srenew(cppopts
, ++ncppopts
);
341 cppopts
[ncppopts
-1] = gmx_strdup(buf
);
349 srenew(cppopts
, ++ncppopts
);
350 cppopts
[ncppopts
-1] = NULL
;
357 find_gb_bondlength(t_params
*plist
, int ai
, int aj
, real
*length
)
364 for (i
= 0; i
< F_NRE
&& !found
; i
++)
368 for (j
= 0; j
< plist
[i
].nr
; j
++)
370 a1
= plist
[i
].param
[j
].a
[0];
371 a2
= plist
[i
].param
[j
].a
[1];
373 if ( (a1
== ai
&& a2
== aj
) || (a1
== aj
&& a2
== ai
))
375 /* Equilibrium bond distance */
376 *length
= plist
[i
].param
[j
].c
[0];
389 find_gb_anglelength(t_params
*plist
, int ai
, int ak
, real
*length
)
391 int i
, j
, a1
, a2
, a3
;
394 int status
, status1
, status2
;
398 for (i
= 0; i
< F_NRE
&& !found
; i
++)
402 for (j
= 0; j
< plist
[i
].nr
; j
++)
404 a1
= plist
[i
].param
[j
].a
[0];
405 a2
= plist
[i
].param
[j
].a
[1];
406 a3
= plist
[i
].param
[j
].a
[2];
408 /* We dont care what the middle atom is, but use it below */
409 if ( (a1
== ai
&& a3
== ak
) || (a1
== ak
&& a3
== ai
) )
411 /* Equilibrium bond distance */
412 a123
= plist
[i
].param
[j
].c
[0];
413 /* Use middle atom to find reference distances r12 and r23 */
414 status1
= find_gb_bondlength(plist
, a1
, a2
, &r12
);
415 status2
= find_gb_bondlength(plist
, a2
, a3
, &r23
);
417 if (status1
== 0 && status2
== 0)
419 /* cosine theorem to get r13 */
420 *length
= std::sqrt(r12
*r12
+r23
*r23
-(2*r12
*r23
*cos(a123
/RAD2DEG
)));
433 generate_gb_exclusion_interactions(t_molinfo
*mi
, gpp_atomtype_t atype
, t_nextnb
*nnb
)
435 int j
, n
, ai
, aj
, ti
, tj
;
440 real radiusi
, radiusj
;
441 real gb_radiusi
, gb_radiusj
;
442 real param_c2
, param_c4
;
448 for (n
= 1; n
<= nnb
->nrex
; n
++)
463 /* Put all higher-order exclusions into 1,4 list so we dont miss them */
470 for (ai
= 0; ai
< nnb
->nr
; ai
++)
472 ti
= at
->atom
[ai
].type
;
473 radiusi
= get_atomtype_radius(ti
, atype
);
474 gb_radiusi
= get_atomtype_gb_radius(ti
, atype
);
476 for (j
= 0; j
< nnb
->nrexcl
[ai
][n
]; j
++)
478 aj
= nnb
->a
[ai
][n
][j
];
480 /* Only add the interactions once */
483 tj
= at
->atom
[aj
].type
;
484 radiusj
= get_atomtype_radius(tj
, atype
);
485 gb_radiusj
= get_atomtype_gb_radius(tj
, atype
);
487 /* There is an exclusion of type "ftype" between atoms ai and aj */
491 /* Reference distance, not used for 1-4 interactions */
495 if (find_gb_bondlength(plist
, ai
, aj
, &distance
) != 0)
497 gmx_fatal(FARGS
, "Cannot find bond length for atoms %d-%d", ai
, aj
);
501 if (find_gb_anglelength(plist
, ai
, aj
, &distance
) != 0)
503 gmx_fatal(FARGS
, "Cannot find length for atoms %d-%d involved in angle", ai
, aj
);
510 /* Assign GB parameters */
512 param
.c
[0] = radiusi
+radiusj
;
513 /* Reference distance distance */
514 param
.c
[1] = distance
;
515 /* Still parameter */
516 param
.c
[2] = param_c2
;
518 param
.c
[3] = gb_radiusi
+gb_radiusj
;
520 param
.c
[4] = param_c4
;
522 /* Add it to the parameter list */
523 add_param_to_list(&plist
[ftype
], ¶m
);
532 static void make_atoms_sys(int nmolb
, const gmx_molblock_t
*molb
,
533 const t_molinfo
*molinfo
,
537 const t_atoms
*mol_atoms
;
542 for (mb
= 0; mb
< nmolb
; mb
++)
544 mol_atoms
= &molinfo
[molb
[mb
].type
].atoms
;
546 srenew(atoms
->atom
, atoms
->nr
+ molb
[mb
].nmol
*mol_atoms
->nr
);
548 for (m
= 0; m
< molb
[mb
].nmol
; m
++)
550 for (a
= 0; a
< mol_atoms
->nr
; a
++)
552 atoms
->atom
[atoms
->nr
++] = mol_atoms
->atom
[a
];
559 static char **read_topol(const char *infile
, const char *outfile
,
560 const char *define
, const char *include
,
562 gpp_atomtype_t atype
,
565 t_molinfo
**intermolecular_interactions
,
567 int *combination_rule
,
572 gmx_molblock_t
**molblock
,
580 char *pline
= NULL
, **title
= NULL
;
581 char line
[STRLEN
], errbuf
[256], comb_str
[256], nb_str
[256];
583 char *dirstr
, *dummy2
;
584 int nrcopies
, nmol
, nmolb
= 0, nscan
, ncombs
, ncopy
;
585 double fLJ
, fQQ
, fPOW
;
586 gmx_molblock_t
*molb
= NULL
;
587 t_molinfo
*mi0
= NULL
;
590 t_nbparam
**nbparam
, **pair
;
592 real fudgeLJ
= -1; /* Multiplication factor to generate 1-4 from LJ */
593 gmx_bool bReadDefaults
, bReadMolType
, bGenPairs
, bWarn_copy_A_B
;
594 double qt
= 0, qBt
= 0; /* total charge */
595 t_bond_atomtype batype
;
597 int dcatt
= -1, nmol_couple
;
598 /* File handling variables */
601 char *tmp_line
= NULL
;
602 char warn_buf
[STRLEN
];
603 const char *floating_point_arithmetic_tip
=
604 "Total charge should normally be an integer. See\n"
605 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
606 "for discussion on how close it should be to an integer.\n";
607 /* We need to open the output file before opening the input file,
608 * because cpp_open_file can change the current working directory.
612 out
= gmx_fio_fopen(outfile
, "w");
619 /* open input file */
620 status
= cpp_open_file(infile
, &handle
, cpp_opts(define
, include
, wi
));
623 gmx_fatal(FARGS
, cpp_error(&handle
, status
));
626 /* some local variables */
627 DS_Init(&DS
); /* directive stack */
628 nmol
= 0; /* no molecules yet... */
629 d
= d_invalid
; /* first thing should be a directive */
630 nbparam
= NULL
; /* The temporary non-bonded matrix */
631 pair
= NULL
; /* The temporary pair interaction matrix */
632 block2
= NULL
; /* the extra exclusions */
635 *reppow
= 12.0; /* Default value for repulsion power */
637 *intermolecular_interactions
= NULL
;
639 /* Init the number of CMAP torsion angles and grid spacing */
640 plist
[F_CMAP
].grid_spacing
= 0;
641 plist
[F_CMAP
].nc
= 0;
643 bWarn_copy_A_B
= bFEP
;
645 batype
= init_bond_atomtype();
646 /* parse the actual file */
647 bReadDefaults
= FALSE
;
649 bReadMolType
= FALSE
;
654 status
= cpp_read_line(&handle
, STRLEN
, line
);
655 done
= (status
== eCPP_EOF
);
658 if (status
!= eCPP_OK
)
660 gmx_fatal(FARGS
, cpp_error(&handle
, status
));
664 fprintf(out
, "%s\n", line
);
667 set_warning_line(wi
, cpp_cur_file(&handle
), cpp_cur_linenr(&handle
));
669 pline
= gmx_strdup(line
);
671 /* Strip trailing '\' from pline, if it exists */
673 if ((sl
> 0) && (pline
[sl
-1] == CONTINUE
))
678 /* build one long line from several fragments - necessary for CMAP */
679 while (continuing(line
))
681 status
= cpp_read_line(&handle
, STRLEN
, line
);
682 set_warning_line(wi
, cpp_cur_file(&handle
), cpp_cur_linenr(&handle
));
684 /* Since we depend on the '\' being present to continue to read, we copy line
685 * to a tmp string, strip the '\' from that string, and cat it to pline
687 tmp_line
= gmx_strdup(line
);
689 sl
= strlen(tmp_line
);
690 if ((sl
> 0) && (tmp_line
[sl
-1] == CONTINUE
))
692 tmp_line
[sl
-1] = ' ';
695 done
= (status
== eCPP_EOF
);
698 if (status
!= eCPP_OK
)
700 gmx_fatal(FARGS
, cpp_error(&handle
, status
));
704 fprintf(out
, "%s\n", line
);
708 srenew(pline
, strlen(pline
)+strlen(tmp_line
)+1);
709 strcat(pline
, tmp_line
);
713 /* skip trailing and leading spaces and comment text */
714 strip_comment (pline
);
717 /* if there is something left... */
718 if ((int)strlen(pline
) > 0)
720 if (pline
[0] == OPENDIR
)
722 /* A directive on this line: copy the directive
723 * without the brackets into dirstr, then
724 * skip spaces and tabs on either side of directive
726 dirstr
= gmx_strdup((pline
+1));
727 if ((dummy2
= strchr (dirstr
, CLOSEDIR
)) != NULL
)
733 if ((newd
= str2dir(dirstr
)) == d_invalid
)
735 sprintf(errbuf
, "Invalid directive %s", dirstr
);
736 warning_error(wi
, errbuf
);
740 /* Directive found */
743 fprintf(debug
, "found directive '%s'\n", dir2str(newd
));
745 if (DS_Check_Order (DS
, newd
))
752 /* we should print here which directives should have
753 been present, and which actually are */
754 gmx_fatal(FARGS
, "%s\nInvalid order for directive %s",
755 cpp_error(&handle
, eCPP_SYNTAX
), dir2str(newd
));
759 if (d
== d_intermolecular_interactions
)
761 if (*intermolecular_interactions
== NULL
)
763 /* We (mis)use the moleculetype processing
764 * to process the intermolecular interactions
765 * by making a "molecule" of the size of the system.
767 snew(*intermolecular_interactions
, 1);
768 init_molinfo(*intermolecular_interactions
);
769 mi0
= *intermolecular_interactions
;
770 make_atoms_sys(nmolb
, molb
, *molinfo
,
777 else if (d
!= d_invalid
)
779 /* Not a directive, just a plain string
780 * use a gigantic switch to decode,
781 * if there is a valid directive!
788 gmx_fatal(FARGS
, "%s\nFound a second defaults directive.\n",
789 cpp_error(&handle
, eCPP_SYNTAX
));
791 bReadDefaults
= TRUE
;
792 nscan
= sscanf(pline
, "%s%s%s%lf%lf%lf",
793 nb_str
, comb_str
, genpairs
, &fLJ
, &fQQ
, &fPOW
);
804 get_nbparm(nb_str
, comb_str
, &nb_funct
, combination_rule
, wi
);
807 bGenPairs
= (gmx_strncasecmp(genpairs
, "Y", 1) == 0);
808 if (nb_funct
!= eNBF_LJ
&& bGenPairs
)
810 gmx_fatal(FARGS
, "Generating pair parameters is only supported with LJ non-bonded interactions");
826 nb_funct
= ifunc_index(d_nonbond_params
, nb_funct
);
830 push_at(symtab
, atype
, batype
, pline
, nb_funct
,
831 &nbparam
, bGenPairs
? &pair
: NULL
, wi
);
835 push_bt(d
, plist
, 2, NULL
, batype
, pline
, wi
);
837 case d_constrainttypes
:
838 push_bt(d
, plist
, 2, NULL
, batype
, pline
, wi
);
843 push_nbt(d
, pair
, atype
, pline
, F_LJ14
, wi
);
847 push_bt(d
, plist
, 2, atype
, NULL
, pline
, wi
);
851 push_bt(d
, plist
, 3, NULL
, batype
, pline
, wi
);
853 case d_dihedraltypes
:
854 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
855 push_dihedraltype(d
, plist
, batype
, pline
, wi
);
858 case d_nonbond_params
:
859 push_nbt(d
, nbparam
, atype
, pline
, nb_funct
, wi
);
864 srenew(block,nblock);
865 srenew(blockinfo,nblock);
866 blk0=&(block[nblock-1]);
867 bi0=&(blockinfo[nblock-1]);
870 push_molt(symtab,bi0,pline);
874 case d_implicit_genborn_params
:
875 push_gb_params(atype
, pline
, wi
);
878 case d_implicit_surface_params
:
879 gmx_fatal(FARGS
, "Implicit surface directive not supported yet.");
883 push_cmaptype(d
, plist
, 5, atype
, batype
, pline
, wi
);
891 if (opts
->couple_moltype
!= NULL
&&
892 (opts
->couple_lam0
== ecouplamNONE
||
893 opts
->couple_lam0
== ecouplamQ
||
894 opts
->couple_lam1
== ecouplamNONE
||
895 opts
->couple_lam1
== ecouplamQ
))
897 dcatt
= add_atomtype_decoupled(symtab
, atype
,
898 &nbparam
, bGenPairs
? &pair
: NULL
);
900 ntype
= get_atomtype_ntypes(atype
);
901 ncombs
= (ntype
*(ntype
+1))/2;
902 generate_nbparams(*combination_rule
, nb_funct
, &(plist
[nb_funct
]), atype
, wi
);
903 ncopy
= copy_nbparams(nbparam
, nb_funct
, &(plist
[nb_funct
]),
905 fprintf(stderr
, "Generated %d of the %d non-bonded parameter combinations\n", ncombs
-ncopy
, ncombs
);
906 free_nbparam(nbparam
, ntype
);
909 gen_pairs(&(plist
[nb_funct
]), &(plist
[F_LJ14
]), fudgeLJ
, *combination_rule
);
910 ncopy
= copy_nbparams(pair
, nb_funct
, &(plist
[F_LJ14
]),
912 fprintf(stderr
, "Generated %d of the %d 1-4 parameter combinations\n", ncombs
-ncopy
, ncombs
);
913 free_nbparam(pair
, ntype
);
915 /* Copy GBSA parameters to atomtype array? */
920 push_molt(symtab
, &nmol
, molinfo
, pline
, wi
);
921 srenew(block2
, nmol
);
922 block2
[nmol
-1].nr
= 0;
923 mi0
= &((*molinfo
)[nmol
-1]);
927 push_atom(symtab
, &(mi0
->cgs
), &(mi0
->atoms
), atype
, pline
, &lastcg
, wi
);
931 push_bond(d
, plist
, mi0
->plist
, &(mi0
->atoms
), atype
, pline
, FALSE
,
932 bGenPairs
, *fudgeQQ
, bZero
, &bWarn_copy_A_B
, wi
);
935 push_bond(d
, plist
, mi0
->plist
, &(mi0
->atoms
), atype
, pline
, FALSE
,
936 FALSE
, 1.0, bZero
, &bWarn_copy_A_B
, wi
);
946 case d_position_restraints
:
947 case d_angle_restraints
:
948 case d_angle_restraints_z
:
949 case d_distance_restraints
:
950 case d_orientation_restraints
:
951 case d_dihedral_restraints
:
954 case d_water_polarization
:
955 case d_thole_polarization
:
956 push_bond(d
, plist
, mi0
->plist
, &(mi0
->atoms
), atype
, pline
, TRUE
,
957 bGenPairs
, *fudgeQQ
, bZero
, &bWarn_copy_A_B
, wi
);
960 push_cmap(d
, plist
, mi0
->plist
, &(mi0
->atoms
), atype
, pline
, wi
);
964 push_vsitesn(d
, mi0
->plist
, &(mi0
->atoms
), pline
, wi
);
967 GMX_ASSERT(block2
, "block2 must always be allocated so exclusions can be processed");
968 if (!block2
[nmol
-1].nr
)
970 init_block2(&(block2
[nmol
-1]), mi0
->atoms
.nr
);
972 push_excl(pline
, &(block2
[nmol
-1]));
976 title
= put_symtab(symtab
, pline
);
983 push_mol(nmol
, *molinfo
, pline
, &whichmol
, &nrcopies
, wi
);
984 mi0
= &((*molinfo
)[whichmol
]);
985 srenew(molb
, nmolb
+1);
986 molb
[nmolb
].type
= whichmol
;
987 molb
[nmolb
].nmol
= nrcopies
;
990 bCouple
= (opts
->couple_moltype
!= NULL
&&
991 (gmx_strcasecmp("system", opts
->couple_moltype
) == 0 ||
992 gmx_strcasecmp(*(mi0
->name
), opts
->couple_moltype
) == 0));
995 nmol_couple
+= nrcopies
;
998 if (mi0
->atoms
.nr
== 0)
1000 gmx_fatal(FARGS
, "Molecule type '%s' contains no atoms",
1004 "Excluding %d bonded neighbours molecule type '%s'\n",
1005 mi0
->nrexcl
, *mi0
->name
);
1006 sum_q(&mi0
->atoms
, nrcopies
, &qt
, &qBt
);
1007 if (!mi0
->bProcessed
)
1010 generate_excl(mi0
->nrexcl
,
1015 merge_excl(&(mi0
->excls
), &(block2
[whichmol
]));
1016 done_block2(&(block2
[whichmol
]));
1017 make_shake(mi0
->plist
, &mi0
->atoms
, opts
->nshake
);
1021 /* nnb contains information about first,2nd,3rd bonded neighbors.
1022 * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
1024 if (bGenborn
== TRUE
)
1026 generate_gb_exclusion_interactions(mi0
, atype
, &nnb
);
1033 convert_moltype_couple(mi0
, dcatt
, *fudgeQQ
,
1034 opts
->couple_lam0
, opts
->couple_lam1
,
1036 nb_funct
, &(plist
[nb_funct
]));
1038 stupid_fill_block(&mi0
->mols
, mi0
->atoms
.nr
, TRUE
);
1039 mi0
->bProcessed
= TRUE
;
1044 fprintf (stderr
, "case: %d\n", (int)d
);
1045 gmx_incons("unknown directive");
1054 status
= cpp_close_file(&handle
);
1055 if (status
!= eCPP_OK
)
1057 gmx_fatal(FARGS
, cpp_error(&handle
, status
));
1062 gmx_fio_fclose(out
);
1065 if (opts
->couple_moltype
)
1067 if (nmol_couple
== 0)
1069 gmx_fatal(FARGS
, "Did not find any molecules of type '%s' for coupling",
1070 opts
->couple_moltype
);
1072 fprintf(stderr
, "Coupling %d copies of molecule type '%s'\n",
1073 nmol_couple
, opts
->couple_moltype
);
1076 /* this is not very clean, but fixes core dump on empty system name */
1079 title
= put_symtab(symtab
, "");
1081 if (fabs(qt
) > 1e-4)
1083 sprintf(warn_buf
, "System has non-zero total charge: %.6f\n%s\n", qt
, floating_point_arithmetic_tip
);
1084 warning_note(wi
, warn_buf
);
1086 if (fabs(qBt
) > 1e-4 && !gmx_within_tol(qBt
, qt
, 1e-6))
1088 sprintf(warn_buf
, "State B has non-zero total charge: %.6f\n%s\n", qBt
, floating_point_arithmetic_tip
);
1089 warning_note(wi
, warn_buf
);
1092 for (i
= 0; i
< nmol
; i
++)
1094 done_block2(&(block2
[i
]));
1098 done_bond_atomtype(&batype
);
1100 if (*intermolecular_interactions
!= NULL
)
1102 sfree(mi0
->atoms
.atom
);
1113 char **do_top(gmx_bool bVerbose
,
1114 const char *topfile
,
1115 const char *topppfile
,
1120 int *combination_rule
,
1121 double *repulsion_power
,
1123 gpp_atomtype_t atype
,
1125 t_molinfo
**molinfo
,
1126 t_molinfo
**intermolecular_interactions
,
1129 gmx_molblock_t
**molblock
,
1133 /* Tmpfile might contain a long path */
1134 const char *tmpfile
;
1139 tmpfile
= topppfile
;
1148 printf("processing topology...\n");
1150 title
= read_topol(topfile
, tmpfile
, opts
->define
, opts
->include
,
1152 nrmols
, molinfo
, intermolecular_interactions
,
1153 plist
, combination_rule
, repulsion_power
,
1154 opts
, fudgeQQ
, nmolblock
, molblock
,
1155 ir
->efep
!= efepNO
, bGenborn
, bZero
, wi
);
1156 if ((*combination_rule
!= eCOMB_GEOMETRIC
) &&
1157 (ir
->vdwtype
== evdwUSER
))
1159 warning(wi
, "Using sigma/epsilon based combination rules with"
1160 " user supplied potential function may produce unwanted"
1168 static void generate_qmexcl_moltype(gmx_moltype_t
*molt
, unsigned char *grpnr
,
1171 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1173 /* generates the exclusions between the individual QM atoms, as
1174 * these interactions should be handled by the QM subroutines and
1175 * not by the gromacs routines
1178 i
, j
, l
, k
= 0, jmax
, qm_max
= 0, qm_nr
= 0, nratoms
= 0, link_nr
= 0, link_max
= 0;
1180 *qm_arr
= NULL
, *link_arr
= NULL
, a1
, a2
, a3
, a4
, ftype
= 0;
1186 *bQMMM
, *blink
, bexcl
;
1188 /* First we search and select the QM atoms in an qm_arr array that
1189 * we use to create the exclusions.
1191 * we take the possibility into account that a user has defined more
1192 * than one QM group:
1194 * for that we also need to do this an ugly work-about just in case
1195 * the QM group contains the entire system...
1197 jmax
= ir
->opts
.ngQM
;
1199 /* we first search for all the QM atoms and put them in an array
1201 for (j
= 0; j
< jmax
; j
++)
1203 for (i
= 0; i
< molt
->atoms
.nr
; i
++)
1205 if (qm_nr
>= qm_max
)
1208 srenew(qm_arr
, qm_max
);
1210 if ((grpnr
? grpnr
[i
] : 0) == j
)
1212 qm_arr
[qm_nr
++] = i
;
1216 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1217 * QM/not QM. We first set all elements to false. Afterwards we use
1218 * the qm_arr to change the elements corresponding to the QM atoms
1221 snew(bQMMM
, molt
->atoms
.nr
);
1222 for (i
= 0; i
< molt
->atoms
.nr
; i
++)
1226 for (i
= 0; i
< qm_nr
; i
++)
1228 bQMMM
[qm_arr
[i
]] = TRUE
;
1231 /* We remove all bonded interactions (i.e. bonds,
1232 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1233 * are removed is as follows: if the interaction invloves 2 atoms,
1234 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1235 * it is removed if at least two of the atoms are QM atoms, if the
1236 * interaction involves 4 atoms, it is removed if there are at least
1237 * 2 QM atoms. Since this routine is called once before any forces
1238 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1239 * be rewritten at this poitn without any problem. 25-9-2002 */
1241 /* first check weter we already have CONNBONDS: */
1242 if (molt
->ilist
[F_CONNBONDS
].nr
!= 0)
1244 fprintf(stderr
, "nr. of CONNBONDS present already: %d\n",
1245 molt
->ilist
[F_CONNBONDS
].nr
/3);
1246 ftype
= molt
->ilist
[F_CONNBONDS
].iatoms
[0];
1247 k
= molt
->ilist
[F_CONNBONDS
].nr
;
1249 /* now we delete all bonded interactions, except the ones describing
1250 * a chemical bond. These are converted to CONNBONDS
1252 for (i
= 0; i
< F_LJ
; i
++)
1254 if (i
== F_CONNBONDS
)
1258 nratoms
= interaction_function
[i
].nratoms
;
1260 while (j
< molt
->ilist
[i
].nr
)
1265 a1
= molt
->ilist
[i
].iatoms
[j
+1];
1266 a2
= molt
->ilist
[i
].iatoms
[j
+2];
1267 bexcl
= (bQMMM
[a1
] && bQMMM
[a2
]);
1268 /* a bonded beteen two QM atoms will be copied to the
1269 * CONNBONDS list, for reasons mentioned above
1271 if (bexcl
&& i
< F_ANGLES
)
1273 srenew(molt
->ilist
[F_CONNBONDS
].iatoms
, k
+3);
1274 molt
->ilist
[F_CONNBONDS
].nr
+= 3;
1275 molt
->ilist
[F_CONNBONDS
].iatoms
[k
++] = ftype
;
1276 molt
->ilist
[F_CONNBONDS
].iatoms
[k
++] = a1
;
1277 molt
->ilist
[F_CONNBONDS
].iatoms
[k
++] = a2
;
1281 a1
= molt
->ilist
[i
].iatoms
[j
+1];
1282 a2
= molt
->ilist
[i
].iatoms
[j
+2];
1283 a3
= molt
->ilist
[i
].iatoms
[j
+3];
1284 bexcl
= ((bQMMM
[a1
] && bQMMM
[a2
]) ||
1285 (bQMMM
[a1
] && bQMMM
[a3
]) ||
1286 (bQMMM
[a2
] && bQMMM
[a3
]));
1289 a1
= molt
->ilist
[i
].iatoms
[j
+1];
1290 a2
= molt
->ilist
[i
].iatoms
[j
+2];
1291 a3
= molt
->ilist
[i
].iatoms
[j
+3];
1292 a4
= molt
->ilist
[i
].iatoms
[j
+4];
1293 bexcl
= ((bQMMM
[a1
] && bQMMM
[a2
] && bQMMM
[a3
]) ||
1294 (bQMMM
[a1
] && bQMMM
[a2
] && bQMMM
[a4
]) ||
1295 (bQMMM
[a1
] && bQMMM
[a3
] && bQMMM
[a4
]) ||
1296 (bQMMM
[a2
] && bQMMM
[a3
] && bQMMM
[a4
]));
1299 gmx_fatal(FARGS
, "no such bonded interactions with %d atoms\n", nratoms
);
1304 /* since the interaction involves QM atoms, these should be
1305 * removed from the MM ilist
1307 molt
->ilist
[i
].nr
-= (nratoms
+1);
1308 for (l
= j
; l
< molt
->ilist
[i
].nr
; l
++)
1310 molt
->ilist
[i
].iatoms
[l
] = molt
->ilist
[i
].iatoms
[l
+(nratoms
+1)];
1315 j
+= nratoms
+1; /* the +1 is for the functype */
1319 /* Now, we search for atoms bonded to a QM atom because we also want
1320 * to exclude their nonbonded interactions with the QM atoms. The
1321 * reason for this is that this interaction is accounted for in the
1322 * linkatoms interaction with the QMatoms and would be counted
1325 for (i
= 0; i
< F_NRE
; i
++)
1330 while (j
< molt
->ilist
[i
].nr
)
1332 a1
= molt
->ilist
[i
].iatoms
[j
+1];
1333 a2
= molt
->ilist
[i
].iatoms
[j
+2];
1334 if ((bQMMM
[a1
] && !bQMMM
[a2
]) || (!bQMMM
[a1
] && bQMMM
[a2
]))
1336 if (link_nr
>= link_max
)
1339 srenew(link_arr
, link_max
);
1343 link_arr
[link_nr
++] = a2
;
1347 link_arr
[link_nr
++] = a1
;
1354 snew(blink
, molt
->atoms
.nr
);
1355 for (i
= 0; i
< molt
->atoms
.nr
; i
++)
1359 for (i
= 0; i
< link_nr
; i
++)
1361 blink
[link_arr
[i
]] = TRUE
;
1363 /* creating the exclusion block for the QM atoms. Each QM atom has
1364 * as excluded elements all the other QMatoms (and itself).
1366 qmexcl
.nr
= molt
->atoms
.nr
;
1367 qmexcl
.nra
= qm_nr
*(qm_nr
+link_nr
)+link_nr
*qm_nr
;
1368 snew(qmexcl
.index
, qmexcl
.nr
+1);
1369 snew(qmexcl
.a
, qmexcl
.nra
);
1371 for (i
= 0; i
< qmexcl
.nr
; i
++)
1373 qmexcl
.index
[i
] = j
;
1376 for (k
= 0; k
< qm_nr
; k
++)
1378 qmexcl
.a
[k
+j
] = qm_arr
[k
];
1380 for (k
= 0; k
< link_nr
; k
++)
1382 qmexcl
.a
[qm_nr
+k
+j
] = link_arr
[k
];
1384 j
+= (qm_nr
+link_nr
);
1388 for (k
= 0; k
< qm_nr
; k
++)
1390 qmexcl
.a
[k
+j
] = qm_arr
[k
];
1395 qmexcl
.index
[qmexcl
.nr
] = j
;
1397 /* and merging with the exclusions already present in sys.
1400 init_block2(&qmexcl2
, molt
->atoms
.nr
);
1401 b_to_b2(&qmexcl
, &qmexcl2
);
1402 merge_excl(&(molt
->excls
), &qmexcl2
);
1403 done_block2(&qmexcl2
);
1405 /* Finally, we also need to get rid of the pair interactions of the
1406 * classical atom bonded to the boundary QM atoms with the QMatoms,
1407 * as this interaction is already accounted for by the QM, so also
1408 * here we run the risk of double counting! We proceed in a similar
1409 * way as we did above for the other bonded interactions: */
1410 for (i
= F_LJ14
; i
< F_COUL14
; i
++)
1412 nratoms
= interaction_function
[i
].nratoms
;
1414 while (j
< molt
->ilist
[i
].nr
)
1416 a1
= molt
->ilist
[i
].iatoms
[j
+1];
1417 a2
= molt
->ilist
[i
].iatoms
[j
+2];
1418 bexcl
= ((bQMMM
[a1
] && bQMMM
[a2
]) ||
1419 (blink
[a1
] && bQMMM
[a2
]) ||
1420 (bQMMM
[a1
] && blink
[a2
]));
1423 /* since the interaction involves QM atoms, these should be
1424 * removed from the MM ilist
1426 molt
->ilist
[i
].nr
-= (nratoms
+1);
1427 for (k
= j
; k
< molt
->ilist
[i
].nr
; k
++)
1429 molt
->ilist
[i
].iatoms
[k
] = molt
->ilist
[i
].iatoms
[k
+(nratoms
+1)];
1434 j
+= nratoms
+1; /* the +1 is for the functype */
1443 } /* generate_qmexcl */
1445 void generate_qmexcl(gmx_mtop_t
*sys
, t_inputrec
*ir
, warninp_t wi
)
1447 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1450 unsigned char *grpnr
;
1451 int mb
, mol
, nat_mol
, i
, nr_mol_with_qm_atoms
= 0;
1452 gmx_molblock_t
*molb
;
1455 grpnr
= sys
->groups
.grpnr
[egcQMMM
];
1457 for (mb
= 0; mb
< sys
->nmolblock
; mb
++)
1459 molb
= &sys
->molblock
[mb
];
1460 nat_mol
= sys
->moltype
[molb
->type
].atoms
.nr
;
1461 for (mol
= 0; mol
< molb
->nmol
; mol
++)
1464 for (i
= 0; i
< nat_mol
; i
++)
1466 if ((grpnr
? grpnr
[i
] : 0) < ir
->opts
.ngQM
)
1473 nr_mol_with_qm_atoms
++;
1476 /* We need to split this molblock */
1479 /* Split the molblock at this molecule */
1481 srenew(sys
->molblock
, sys
->nmolblock
);
1482 for (i
= sys
->nmolblock
-2; i
>= mb
; i
--)
1484 sys
->molblock
[i
+1] = sys
->molblock
[i
];
1486 sys
->molblock
[mb
].nmol
= mol
;
1487 sys
->molblock
[mb
+1].nmol
-= mol
;
1489 molb
= &sys
->molblock
[mb
];
1493 /* Split the molblock after this molecule */
1495 srenew(sys
->molblock
, sys
->nmolblock
);
1496 molb
= &sys
->molblock
[mb
];
1497 for (i
= sys
->nmolblock
-2; i
>= mb
; i
--)
1499 sys
->molblock
[i
+1] = sys
->molblock
[i
];
1501 sys
->molblock
[mb
].nmol
= 1;
1502 sys
->molblock
[mb
+1].nmol
-= 1;
1505 /* Add a moltype for the QMMM molecule */
1507 srenew(sys
->moltype
, sys
->nmoltype
);
1508 /* Copy the moltype struct */
1509 sys
->moltype
[sys
->nmoltype
-1] = sys
->moltype
[molb
->type
];
1510 /* Copy the exclusions to a new array, since this is the only
1511 * thing that needs to be modified for QMMM.
1513 copy_blocka(&sys
->moltype
[molb
->type
].excls
,
1514 &sys
->moltype
[sys
->nmoltype
-1].excls
);
1515 /* Set the molecule type for the QMMM molblock */
1516 molb
->type
= sys
->nmoltype
- 1;
1518 generate_qmexcl_moltype(&sys
->moltype
[molb
->type
], grpnr
, ir
);
1526 if (nr_mol_with_qm_atoms
> 1)
1528 /* generate a warning is there are QM atoms in different
1529 * topologies. In this case it is not possible at this stage to
1530 * mutualy exclude the non-bonded interactions via the
1531 * exclusions (AFAIK). Instead, the user is advised to use the
1532 * energy group exclusions in the mdp file
1535 "\nThe QM subsystem is divided over multiple topologies. "
1536 "The mutual non-bonded interactions cannot be excluded. "
1537 "There are two ways to achieve this:\n\n"
1538 "1) merge the topologies, such that the atoms of the QM "
1539 "subsystem are all present in one single topology file. "
1540 "In this case this warning will dissappear\n\n"
1541 "2) exclude the non-bonded interactions explicitly via the "
1542 "energygrp-excl option in the mdp file. if this is the case "
1543 "this warning may be ignored"