1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
41 #include <sys/types.h>
60 #include "gmx_fatal.h"
62 #include "vsite_parm.h"
68 #include "gpp_nextnb.h"
72 #include "gpp_bond_atomtype.h"
76 #define CPPMARK '#' /* mark from cpp */
77 #define OPENDIR '[' /* starting sign for directive */
78 #define CLOSEDIR ']' /* ending sign for directive */
81 printf("line: %d, maxavail: %d\n",__LINE__,maxavail()); \
84 static void free_nbparam(t_nbparam **param,int nr)
93 static int copy_nbparams(t_nbparam **param,int ftype,t_params *plist,int nr)
103 if (param[i][j].bSet) {
104 for(f=0; f<nrfp; f++) {
105 plist->param[nr*i+j].c[f] = param[i][j].c[f];
106 plist->param[nr*j+i].c[f] = param[i][j].c[f];
114 static void gen_pairs(t_params *nbs,t_params *pairs,real fudge, int comb, gmx_bool bVerbose)
116 int i,j,ntp,nrfp,nrfpA,nrfpB,nnn;
121 nrfpA = interaction_function[F_LJ14].nrfpA;
122 nrfpB = interaction_function[F_LJ14].nrfpB;
125 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
126 gmx_incons("Number of force parameters in gen_pairs wrong");
128 fprintf(stderr,"Generating 1-4 interactions: fudge = %g\n",fudge);
130 fprintf(debug,"Fudge factor for 1-4 interactions: %g\n",fudge);
131 fprintf(debug,"Holy Cow! there are %d types\n",ntp);
133 snew(pairs->param,pairs->nr);
134 for(i=0; (i<ntp); i++) {
136 pairs->param[i].a[0] = i / nnn;
137 pairs->param[i].a[1] = i % nnn;
138 /* Copy normal and FEP parameters and multiply by fudge factor */
142 for(j=0; (j<nrfp); j++) {
143 /* If we are using sigma/epsilon values, only the epsilon values
144 * should be scaled, but not sigma.
145 * The sigma values have even indices 0,2, etc.
147 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2==0))
152 pairs->param[i].c[j]=scaling*nbs->param[i].c[j];
153 pairs->param[i].c[nrfp+j]=scaling*nbs->param[i].c[j];
158 double check_mol(gmx_mtop_t *mtop,warninp_t wi)
166 /* Check mass and charge */
169 for(mb=0; mb<mtop->nmoltype; mb++) {
170 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
171 nmol = mtop->molblock[mb].nmol;
172 for (i=0; (i<atoms->nr); i++) {
173 q += nmol*atoms->atom[i].q;
174 m = atoms->atom[i].m;
175 pt = atoms->atom[i].ptype;
176 /* If the particle is an atom or a nucleus it must have a mass,
177 * else, if it is a shell, a vsite or a bondshell it can have mass zero
179 if ((m <= 0.0) && ((pt == eptAtom) || (pt == eptNucleus))) {
180 ri = atoms->atom[i].resind;
181 sprintf(buf,"atom %s (Res %s-%d) has mass %g\n",
182 *(atoms->atomname[i]),
183 *(atoms->resinfo[ri].name),
184 atoms->resinfo[ri].nr,
186 warning_error(wi,buf);
188 if ((m!=0) && (pt == eptVSite)) {
189 ri = atoms->atom[i].resind;
190 sprintf(buf,"virtual site %s (Res %s-%d) has non-zero mass %g\n"
191 " Check your topology.\n",
192 *(atoms->atomname[i]),
193 *(atoms->resinfo[ri].name),
194 atoms->resinfo[ri].nr,
196 warning_error(wi,buf);
197 /* The following statements make LINCS break! */
198 /* atoms->atom[i].m=0; */
205 static void sum_q(t_atoms *atoms,int n,double *qt,double *qBt)
213 for (i=0; i<atoms->nr; i++)
215 qmolA += atoms->atom[i].q;
216 qmolB += atoms->atom[i].qB;
218 /* Unfortunately an absolute comparison,
219 * but this avoids unnecessary warnings and gmx-users mails.
221 if (fabs(qmolA) >= 1e-6 || fabs(qmolB) >= 1e-6)
228 static void get_nbparm(char *nb_str,char *comb_str,int *nb,int *comb,
232 char warn_buf[STRLEN];
235 for(i=1; (i<eNBF_NR); i++)
236 if (gmx_strcasecmp(nb_str,enbf_names[i]) == 0)
239 *nb = strtol(nb_str,NULL,10);
240 if ((*nb < 1) || (*nb >= eNBF_NR)) {
241 sprintf(warn_buf,"Invalid nonbond function selector '%s' using %s",
242 nb_str,enbf_names[1]);
243 warning_error(wi,warn_buf);
247 for(i=1; (i<eCOMB_NR); i++)
248 if (gmx_strcasecmp(comb_str,ecomb_names[i]) == 0)
251 *comb = strtol(comb_str,NULL,10);
252 if ((*comb < 1) || (*comb >= eCOMB_NR)) {
253 sprintf(warn_buf,"Invalid combination rule selector '%s' using %s",
254 comb_str,ecomb_names[1]);
255 warning_error(wi,warn_buf);
260 static char ** cpp_opts(const char *define,const char *include,
266 const char *cppadds[2];
267 char **cppopts = NULL;
268 const char *option[2] = { "-D","-I" };
269 const char *nopt[2] = { "define", "include" };
273 char warn_buf[STRLEN];
276 cppadds[1] = include;
277 for(n=0; (n<2); n++) {
280 while(*ptr != '\0') {
281 while((*ptr != '\0') && isspace(*ptr))
284 while((*rptr != '\0') && !isspace(*rptr))
289 strncpy(buf,ptr,len);
290 if (strstr(ptr,option[n]) != ptr) {
291 set_warning_line(wi,"mdp file",-1);
292 sprintf(warn_buf,"Malformed %s option %s",nopt[n],buf);
293 warning(wi,warn_buf);
296 srenew(cppopts,++ncppopts);
297 cppopts[ncppopts-1] = strdup(buf);
305 srenew(cppopts,++ncppopts);
306 cppopts[ncppopts-1] = NULL;
313 find_gb_bondlength(t_params *plist,int ai,int aj, real *length)
320 for(i=0;i<F_NRE && !found;i++)
324 for(j=0;j<plist[i].nr; j++)
326 a1=plist[i].param[j].a[0];
327 a2=plist[i].param[j].a[1];
329 if( (a1==ai && a2==aj) || (a1==aj && a2==ai))
331 /* Equilibrium bond distance */
332 *length=plist[i].param[j].c[0];
345 find_gb_anglelength(t_params *plist,int ai,int ak, real *length)
350 int status,status1,status2;
354 for(i=0;i<F_NRE && !found;i++)
358 for(j=0;j<plist[i].nr; j++)
360 a1=plist[i].param[j].a[0];
361 a2=plist[i].param[j].a[1];
362 a3=plist[i].param[j].a[2];
364 /* We dont care what the middle atom is, but use it below */
365 if( (a1==ai && a3==ak) || (a1==ak && a3==ai) )
367 /* Equilibrium bond distance */
368 a123 = plist[i].param[j].c[0];
369 /* Use middle atom to find reference distances r12 and r23 */
370 status1 = find_gb_bondlength(plist,a1,a2,&r12);
371 status2 = find_gb_bondlength(plist,a2,a3,&r23);
373 if(status1==0 && status2==0)
375 /* cosine theorem to get r13 */
376 *length=sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG)));
389 generate_gb_exclusion_interactions(t_molinfo *mi,gpp_atomtype_t atype,t_nextnb *nnb)
391 int i,j,k,n,ai,aj,ti,tj;
397 real radiusi,radiusj;
398 real gb_radiusi,gb_radiusj;
399 real param_c2,param_c4;
405 for(n=1;n<=nnb->nrex;n++)
420 /* Put all higher-order exclusions into 1,4 list so we dont miss them */
427 for(ai=0;ai<nnb->nr;ai++)
429 ti = at->atom[ai].type;
430 radiusi = get_atomtype_radius(ti,atype);
431 gb_radiusi = get_atomtype_gb_radius(ti,atype);
433 for(j=0;j<nnb->nrexcl[ai][n];j++)
435 aj = nnb->a[ai][n][j];
437 /* Only add the interactions once */
440 tj = at->atom[aj].type;
441 radiusj = get_atomtype_radius(tj,atype);
442 gb_radiusj = get_atomtype_gb_radius(tj,atype);
444 /* There is an exclusion of type "ftype" between atoms ai and aj */
448 /* Reference distance, not used for 1-4 interactions */
452 if(find_gb_bondlength(plist,ai,aj,&distance)!=0)
454 gmx_fatal(FARGS,"Cannot find bond length for atoms %d-%d",ai,aj);
458 if(find_gb_anglelength(plist,ai,aj,&distance)!=0)
460 gmx_fatal(FARGS,"Cannot find length for atoms %d-%d involved in angle",ai,aj);
467 /* Assign GB parameters */
469 param.c[0] = radiusi+radiusj;
470 /* Reference distance distance */
471 param.c[1] = distance;
472 /* Still parameter */
473 param.c[2] = param_c2;
475 param.c[3] = gb_radiusi+gb_radiusj;
477 param.c[4] = param_c4;
479 /* Add it to the parameter list */
480 add_param_to_list(&plist[ftype],¶m);
491 static char **read_topol(const char *infile,const char *outfile,
492 const char *define,const char *include,
494 gpp_atomtype_t atype,
498 int *combination_rule,
503 gmx_molblock_t **molblock,
511 int i,sl,nb_funct,comb;
512 char *pline=NULL,**title=NULL;
513 char line[STRLEN],errbuf[256],comb_str[256],nb_str[256];
515 char *dirstr,*dummy2;
516 int nrcopies,nmol,nmolb=0,nscan,ncombs,ncopy;
518 gmx_molblock_t *molb=NULL;
519 t_topology *block=NULL;
523 t_nbparam **nbparam,**pair;
525 real fudgeLJ=-1; /* Multiplication factor to generate 1-4 from LJ */
526 gmx_bool bReadDefaults,bReadMolType,bGenPairs,bWarn_copy_A_B;
527 double qt=0,qBt=0; /* total charge */
528 t_bond_atomtype batype;
530 int dcatt=-1,nmol_couple;
531 /* File handling variables */
535 char warn_buf[STRLEN];
536 const char *floating_point_arithmetic_tip =
537 "Total charge should normally be an integer. See\n"
538 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
539 "for discussion on how close it should be to an integer.\n";
540 /* We need to open the output file before opening the input file,
541 * because cpp_open_file can change the current working directory.
544 out = gmx_fio_fopen(outfile,"w");
549 /* open input file */
550 status = cpp_open_file(infile,&handle,cpp_opts(define,include,infile,wi));
552 gmx_fatal(FARGS,cpp_error(&handle,status));
554 /* some local variables */
555 DS_Init(&DS); /* directive stack */
556 nmol = 0; /* no molecules yet... */
557 d = d_invalid; /* first thing should be a directive */
558 nbparam = NULL; /* The temporary non-bonded matrix */
559 pair = NULL; /* The temporary pair interaction matrix */
560 block2 = NULL; /* the extra exclusions */
562 *reppow = 12.0; /* Default value for repulsion power */
566 /* Init the number of CMAP torsion angles and grid spacing */
567 plist->grid_spacing = 0;
570 bWarn_copy_A_B = bFEP;
572 batype = init_bond_atomtype();
573 /* parse the actual file */
574 bReadDefaults = FALSE;
576 bReadMolType = FALSE;
580 status = cpp_read_line(&handle,STRLEN,line);
581 done = (status == eCPP_EOF);
583 if (status != eCPP_OK)
584 gmx_fatal(FARGS,cpp_error(&handle,status));
586 fprintf(out,"%s\n",line);
588 set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
590 pline = strdup(line);
592 /* Strip trailing '\' from pline, if it exists */
594 if ((sl > 0) && (pline[sl-1] == CONTINUE)) {
598 /* build one long line from several fragments - necessary for CMAP */
599 while (continuing(line))
601 status = cpp_read_line(&handle,STRLEN,line);
602 set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
604 /* Since we depend on the '\' being present to continue to read, we copy line
605 * to a tmp string, strip the '\' from that string, and cat it to pline
607 tmp_line=strdup(line);
609 sl = strlen(tmp_line);
610 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE)) {
611 tmp_line[sl-1] = ' ';
614 done = (status == eCPP_EOF);
616 if (status != eCPP_OK)
617 gmx_fatal(FARGS,cpp_error(&handle,status));
619 fprintf(out,"%s\n",line);
622 srenew(pline,strlen(pline)+strlen(tmp_line)+1);
623 strcat(pline,tmp_line);
627 /* skip trailing and leading spaces and comment text */
628 strip_comment (pline);
631 /* if there is something left... */
632 if ((int)strlen(pline) > 0) {
633 if (pline[0] == OPENDIR) {
634 /* A directive on this line: copy the directive
635 * without the brackets into dirstr, then
636 * skip spaces and tabs on either side of directive
638 dirstr = strdup((pline+1));
639 if ((dummy2 = strchr (dirstr,CLOSEDIR)) != NULL)
643 if ((newd = str2dir(dirstr)) == d_invalid) {
644 sprintf(errbuf,"Invalid directive %s",dirstr);
645 warning_error(wi,errbuf);
648 /* Directive found */
650 fprintf(debug,"found directive '%s'\n",dir2str(newd));
651 if (DS_Check_Order (DS,newd)) {
656 /* we should print here which directives should have
657 been present, and which actually are */
658 gmx_fatal(FARGS,"%s\nInvalid order for directive %s",
659 cpp_error(&handle,eCPP_SYNTAX),dir2str(newd));
665 else if (d != d_invalid) {
666 /* Not a directive, just a plain string
667 * use a gigantic switch to decode,
668 * if there is a valid directive!
673 gmx_fatal(FARGS,"%s\nFound a second defaults directive.\n",
674 cpp_error(&handle,eCPP_SYNTAX));
675 bReadDefaults = TRUE;
676 nscan = sscanf(pline,"%s%s%s%lf%lf%lf",
677 nb_str,comb_str,genpairs,&fLJ,&fQQ,&fPOW);
685 get_nbparm(nb_str,comb_str,&nb_funct,&comb,wi);
686 *combination_rule = comb;
688 bGenPairs = (gmx_strncasecmp(genpairs,"Y",1) == 0);
689 if (nb_funct != eNBF_LJ && bGenPairs) {
690 gmx_fatal(FARGS,"Generating pair parameters is only supported with LJ non-bonded interactions");
700 nb_funct = ifunc_index(d_nonbond_params,nb_funct);
704 push_at(symtab,atype,batype,pline,nb_funct,
705 &nbparam,bGenPairs ? &pair : NULL,wi);
709 push_bt(d,plist,2,NULL,batype,pline,wi);
711 case d_constrainttypes:
712 push_bt(d,plist,2,NULL,batype,pline,wi);
716 push_nbt(d,pair,atype,pline,F_LJ14,wi);
718 push_bt(d,plist,2,atype,NULL,pline,wi);
721 push_bt(d,plist,3,NULL,batype,pline,wi);
723 case d_dihedraltypes:
724 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
725 push_dihedraltype(d,plist,batype,pline,wi);
728 case d_nonbond_params:
729 push_nbt(d,nbparam,atype,pline,nb_funct,wi);
734 srenew(block,nblock);
735 srenew(blockinfo,nblock);
736 blk0=&(block[nblock-1]);
737 bi0=&(blockinfo[nblock-1]);
740 push_molt(symtab,bi0,pline);
744 case d_implicit_genborn_params:
745 push_gb_params(atype,pline,wi);
748 case d_implicit_surface_params:
749 gmx_fatal(FARGS,"Implicit surface directive not supported yet.");
753 push_cmaptype(d, plist, 5, atype, batype,pline,wi);
756 case d_moleculetype: {
759 if (opts->couple_moltype != NULL &&
760 (opts->couple_lam0 == ecouplamNONE ||
761 opts->couple_lam0 == ecouplamQ ||
762 opts->couple_lam1 == ecouplamNONE ||
763 opts->couple_lam1 == ecouplamQ))
765 dcatt = add_atomtype_decoupled(symtab,atype,
766 &nbparam,bGenPairs?&pair:NULL);
768 ntype = get_atomtype_ntypes(atype);
769 ncombs = (ntype*(ntype+1))/2;
770 generate_nbparams(comb,nb_funct,&(plist[nb_funct]),atype,wi);
771 ncopy = copy_nbparams(nbparam,nb_funct,&(plist[nb_funct]),
773 fprintf(stderr,"Generated %d of the %d non-bonded parameter combinations\n",ncombs-ncopy,ncombs);
774 free_nbparam(nbparam,ntype);
776 gen_pairs(&(plist[nb_funct]),&(plist[F_LJ14]),fudgeLJ,comb,bVerbose);
777 ncopy = copy_nbparams(pair,nb_funct,&(plist[F_LJ14]),
779 fprintf(stderr,"Generated %d of the %d 1-4 parameter combinations\n",ncombs-ncopy,ncombs);
780 free_nbparam(pair,ntype);
782 /* Copy GBSA parameters to atomtype array? */
787 push_molt(symtab,&nmol,molinfo,pline,wi);
790 mi0=&((*molinfo)[nmol-1]);
794 push_atom(symtab,&(mi0->cgs),&(mi0->atoms),atype,pline,&lastcg,wi);
798 push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,FALSE,
799 bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
809 case d_position_restraints:
810 case d_angle_restraints:
811 case d_angle_restraints_z:
812 case d_distance_restraints:
813 case d_orientation_restraints:
814 case d_dihedral_restraints:
817 case d_water_polarization:
818 case d_thole_polarization:
819 push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,TRUE,
820 bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
823 push_cmap(d,plist,mi0->plist,&(mi0->atoms),atype,pline,
828 push_vsitesn(d,plist,mi0->plist,&(mi0->atoms),atype,pline,wi);
832 if (!block2[nmol-1].nr)
833 init_block2(&(block2[nmol-1]),mi0->atoms.nr);
834 push_excl(pline,&(block2[nmol-1]));
838 title=put_symtab(symtab,pline);
844 push_mol(nmol,*molinfo,pline,&whichmol,&nrcopies,wi);
845 mi0=&((*molinfo)[whichmol]);
846 srenew(molb,nmolb+1);
847 molb[nmolb].type = whichmol;
848 molb[nmolb].nmol = nrcopies;
851 bCouple = (opts->couple_moltype != NULL &&
852 (gmx_strcasecmp("system" ,opts->couple_moltype) == 0 ||
853 gmx_strcasecmp(*(mi0->name),opts->couple_moltype) == 0));
855 nmol_couple += nrcopies;
858 if (mi0->atoms.nr == 0) {
859 gmx_fatal(FARGS,"Molecule type '%s' contains no atoms",
863 "Excluding %d bonded neighbours molecule type '%s'\n",
864 mi0->nrexcl,*mi0->name);
865 sum_q(&mi0->atoms,nrcopies,&qt,&qBt);
866 if (!mi0->bProcessed)
869 generate_excl(mi0->nrexcl,
874 merge_excl(&(mi0->excls),&(block2[whichmol]));
875 done_block2(&(block2[whichmol]));
876 make_shake(mi0->plist,&mi0->atoms,atype,opts->nshake);
880 /* nnb contains information about first,2nd,3rd bonded neighbors.
881 * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
885 generate_gb_exclusion_interactions(mi0,atype,&nnb);
892 convert_moltype_couple(mi0,dcatt,*fudgeQQ,
893 opts->couple_lam0,opts->couple_lam1,
895 nb_funct,&(plist[nb_funct]));
897 stupid_fill_block(&mi0->mols,mi0->atoms.nr,TRUE);
898 mi0->bProcessed=TRUE;
903 fprintf (stderr,"case: %d\n",d);
912 status = cpp_close_file(&handle);
913 if (status != eCPP_OK)
914 gmx_fatal(FARGS,cpp_error(&handle,status));
918 if (opts->couple_moltype) {
919 if (nmol_couple == 0) {
920 gmx_fatal(FARGS,"Did not find any molecules of type '%s' for coupling",
921 opts->couple_moltype);
923 fprintf(stderr,"Coupling %d copies of molecule type '%s'\n",
924 nmol_couple,opts->couple_moltype);
927 /* this is not very clean, but fixes core dump on empty system name */
929 title=put_symtab(symtab,"");
930 if (fabs(qt) > 1e-4) {
931 sprintf(warn_buf,"System has non-zero total charge: %.6f\n%s\n",qt,floating_point_arithmetic_tip);
932 warning_note(wi,warn_buf);
934 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt,qt,1e-6)) {
935 sprintf(warn_buf,"State B has non-zero total charge: %.6f\n%s\n",qBt,floating_point_arithmetic_tip);
936 warning_note(wi,warn_buf);
939 for(i=0; i<nmol; i++)
940 done_block2(&(block2[i]));
943 done_bond_atomtype(&batype);
953 char **do_top(gmx_bool bVerbose,
955 const char *topppfile,
960 int *combination_rule,
961 double *repulsion_power,
963 gpp_atomtype_t atype,
968 gmx_molblock_t **molblock,
972 /* Tmpfile might contain a long path */
987 printf("processing topology...\n");
989 title = read_topol(topfile,tmpfile,opts->define,opts->include,
990 symtab,atype,nrmols,molinfo,
991 plist,combination_rule,repulsion_power,
992 opts,fudgeQQ,nmolblock,molblock,
993 ir->efep!=efepNO,bGenborn,bZero,bVerbose,
995 if ((*combination_rule != eCOMB_GEOMETRIC) &&
996 (ir->vdwtype == evdwUSER))
998 warning(wi,"Using sigma/epsilon based combination rules with"
999 " user supplied potential function may produce unwanted"
1007 static void generate_qmexcl_moltype(gmx_moltype_t *molt,unsigned char *grpnr,
1010 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1012 /* generates the exclusions between the individual QM atoms, as
1013 * these interactions should be handled by the QM subroutines and
1014 * not by the gromacs routines
1017 i,j,l,k=0,jmax,qm_max=0,qm_nr=0,nratoms=0,link_nr=0,link_max=0;
1019 *qm_arr=NULL,*link_arr=NULL,a1,a2,a3,a4,ftype=0;
1025 *bQMMM,*blink,bexcl;
1027 /* First we search and select the QM atoms in an qm_arr array that
1028 * we use to create the exclusions.
1030 * we take the possibility into account that a user has defined more
1031 * than one QM group:
1033 * for that we also need to do this an ugly work-about just in case
1034 * the QM group contains the entire system...
1036 jmax = ir->opts.ngQM;
1038 /* we first search for all the QM atoms and put them in an array
1040 for(j=0;j<jmax;j++){
1041 for(i=0;i<molt->atoms.nr;i++){
1044 srenew(qm_arr,qm_max);
1046 if ((grpnr ? grpnr[i] : 0) == j){
1047 qm_arr[qm_nr++] = i;
1051 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1052 * QM/not QM. We first set all elements to false. Afterwards we use
1053 * the qm_arr to change the elements corresponding to the QM atoms
1056 snew(bQMMM,molt->atoms.nr);
1057 for (i=0;i<molt->atoms.nr;i++)
1059 for (i=0;i<qm_nr;i++)
1060 bQMMM[qm_arr[i]]=TRUE;
1062 /* We remove all bonded interactions (i.e. bonds,
1063 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1064 * are removed is as follows: if the interaction invloves 2 atoms,
1065 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1066 * it is removed if at least two of the atoms are QM atoms, if the
1067 * interaction involves 4 atoms, it is removed if there are at least
1068 * 2 QM atoms. Since this routine is called once before any forces
1069 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1070 * be rewritten at this poitn without any problem. 25-9-2002 */
1072 /* first check weter we already have CONNBONDS: */
1073 if (molt->ilist[F_CONNBONDS].nr != 0){
1074 fprintf(stderr,"nr. of CONNBONDS present already: %d\n",
1075 molt->ilist[F_CONNBONDS].nr/3);
1076 ftype = molt->ilist[F_CONNBONDS].iatoms[0];
1077 k = molt->ilist[F_CONNBONDS].nr;
1079 /* now we delete all bonded interactions, except the ones describing
1080 * a chemical bond. These are converted to CONNBONDS
1082 for (i=0;i<F_LJ;i++){
1085 nratoms = interaction_function[i].nratoms;
1087 while (j<molt->ilist[i].nr){
1091 a1 = molt->ilist[i].iatoms[j+1];
1092 a2 = molt->ilist[i].iatoms[j+2];
1093 bexcl = (bQMMM[a1] && bQMMM[a2]);
1094 /* a bonded beteen two QM atoms will be copied to the
1095 * CONNBONDS list, for reasons mentioned above
1097 if(bexcl && i<F_ANGLES){
1098 srenew(molt->ilist[F_CONNBONDS].iatoms,k+3);
1099 molt->ilist[F_CONNBONDS].nr += 3;
1100 molt->ilist[F_CONNBONDS].iatoms[k++] = ftype;
1101 molt->ilist[F_CONNBONDS].iatoms[k++] = a1;
1102 molt->ilist[F_CONNBONDS].iatoms[k++] = a2;
1106 a1 = molt->ilist[i].iatoms[j+1];
1107 a2 = molt->ilist[i].iatoms[j+2];
1108 a3 = molt->ilist[i].iatoms[j+3];
1109 bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1110 (bQMMM[a1] && bQMMM[a3]) ||
1111 (bQMMM[a2] && bQMMM[a3]));
1114 a1 = molt->ilist[i].iatoms[j+1];
1115 a2 = molt->ilist[i].iatoms[j+2];
1116 a3 = molt->ilist[i].iatoms[j+3];
1117 a4 = molt->ilist[i].iatoms[j+4];
1118 bexcl = ((bQMMM[a1] && bQMMM[a2] && bQMMM[a3]) ||
1119 (bQMMM[a1] && bQMMM[a2] && bQMMM[a4]) ||
1120 (bQMMM[a1] && bQMMM[a3] && bQMMM[a4]) ||
1121 (bQMMM[a2] && bQMMM[a3] && bQMMM[a4]));
1124 gmx_fatal(FARGS,"no such bonded interactions with %d atoms\n",nratoms);
1127 /* since the interaction involves QM atoms, these should be
1128 * removed from the MM ilist
1130 molt->ilist[i].nr -= (nratoms+1);
1131 for (l=j;l<molt->ilist[i].nr;l++)
1132 molt->ilist[i].iatoms[l] = molt->ilist[i].iatoms[l+(nratoms+1)];
1134 j += nratoms+1; /* the +1 is for the functype */
1138 /* Now, we search for atoms bonded to a QM atom because we also want
1139 * to exclude their nonbonded interactions with the QM atoms. The
1140 * reason for this is that this interaction is accounted for in the
1141 * linkatoms interaction with the QMatoms and would be counted
1144 for(i=0;i<F_NRE;i++){
1147 while(j<molt->ilist[i].nr){
1148 a1 = molt->ilist[i].iatoms[j+1];
1149 a2 = molt->ilist[i].iatoms[j+2];
1150 if((bQMMM[a1] && !bQMMM[a2])||(!bQMMM[a1] && bQMMM[a2])){
1151 if(link_nr>=link_max){
1153 srenew(link_arr,link_max);
1156 link_arr[link_nr++] = a2;
1158 link_arr[link_nr++] = a1;
1165 snew(blink,molt->atoms.nr);
1166 for (i=0;i<molt->atoms.nr;i++)
1168 for (i=0;i<link_nr;i++)
1169 blink[link_arr[i]]=TRUE;
1170 /* creating the exclusion block for the QM atoms. Each QM atom has
1171 * as excluded elements all the other QMatoms (and itself).
1173 qmexcl.nr = molt->atoms.nr;
1174 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1175 snew(qmexcl.index,qmexcl.nr+1);
1176 snew(qmexcl.a,qmexcl.nra);
1178 for(i=0;i<qmexcl.nr;i++){
1181 for(k=0;k<qm_nr;k++){
1182 qmexcl.a[k+j] = qm_arr[k];
1184 for(k=0;k<link_nr;k++){
1185 qmexcl.a[qm_nr+k+j] = link_arr[k];
1190 for(k=0;k<qm_nr;k++){
1191 qmexcl.a[k+j] = qm_arr[k];
1196 qmexcl.index[qmexcl.nr]=j;
1198 /* and merging with the exclusions already present in sys.
1201 init_block2(&qmexcl2,molt->atoms.nr);
1202 b_to_b2(&qmexcl, &qmexcl2);
1203 merge_excl(&(molt->excls),&qmexcl2);
1204 done_block2(&qmexcl2);
1206 /* Finally, we also need to get rid of the pair interactions of the
1207 * classical atom bonded to the boundary QM atoms with the QMatoms,
1208 * as this interaction is already accounted for by the QM, so also
1209 * here we run the risk of double counting! We proceed in a similar
1210 * way as we did above for the other bonded interactions: */
1211 for (i=F_LJ14;i<F_COUL14;i++){
1212 nratoms = interaction_function[i].nratoms;
1214 while (j<molt->ilist[i].nr){
1216 a1 = molt->ilist[i].iatoms[j+1];
1217 a2 = molt->ilist[i].iatoms[j+2];
1218 bexcl = ((bQMMM[a1] && bQMMM[a2])||
1219 (blink[a1] && bQMMM[a2])||
1220 (bQMMM[a1] && blink[a2]));
1222 /* since the interaction involves QM atoms, these should be
1223 * removed from the MM ilist
1225 molt->ilist[i].nr -= (nratoms+1);
1226 for (k=j;k<molt->ilist[i].nr;k++)
1227 molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1229 j += nratoms+1; /* the +1 is for the functype */
1238 } /* generate_qmexcl */
1240 void generate_qmexcl(gmx_mtop_t *sys,t_inputrec *ir,warninp_t wi)
1242 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1245 unsigned char *grpnr;
1246 int mb,mol,nat_mol,i,nr_mol_with_qm_atoms=0;
1247 gmx_molblock_t *molb;
1250 grpnr = sys->groups.grpnr[egcQMMM];
1252 for(mb=0; mb<sys->nmolblock; mb++) {
1253 molb = &sys->molblock[mb];
1254 nat_mol = sys->moltype[molb->type].atoms.nr;
1255 for(mol=0; mol<molb->nmol; mol++) {
1257 for(i=0; i<nat_mol; i++) {
1258 if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM) {
1263 nr_mol_with_qm_atoms++;
1264 if (molb->nmol > 1) {
1265 /* We need to split this molblock */
1267 /* Split the molblock at this molecule */
1269 srenew(sys->molblock,sys->nmolblock);
1270 for(i=sys->nmolblock-2; i >= mb;i--) {
1271 sys->molblock[i+1] = sys->molblock[i];
1273 sys->molblock[mb ].nmol = mol;
1274 sys->molblock[mb+1].nmol -= mol;
1276 molb = &sys->molblock[mb];
1278 if (molb->nmol > 1) {
1279 /* Split the molblock after this molecule */
1281 srenew(sys->molblock,sys->nmolblock);
1282 molb = &sys->molblock[mb];
1283 for(i=sys->nmolblock-2; i >= mb;i--) {
1284 sys->molblock[i+1] = sys->molblock[i];
1286 sys->molblock[mb ].nmol = 1;
1287 sys->molblock[mb+1].nmol -= 1;
1290 /* Add a moltype for the QMMM molecule */
1292 srenew(sys->moltype,sys->nmoltype);
1293 /* Copy the moltype struct */
1294 sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
1295 /* Copy the exclusions to a new array, since this is the only
1296 * thing that needs to be modified for QMMM.
1298 copy_blocka(&sys->moltype[molb->type ].excls,
1299 &sys->moltype[sys->nmoltype-1].excls);
1300 /* Set the molecule type for the QMMM molblock */
1301 molb->type = sys->nmoltype - 1;
1303 generate_qmexcl_moltype(&sys->moltype[molb->type],grpnr,ir);
1310 if(nr_mol_with_qm_atoms>1){
1311 /* generate a warning is there are QM atoms in different
1312 * topologies. In this case it is not possible at this stage to
1313 * mutualy exclude the non-bonded interactions via the
1314 * exclusions (AFAIK). Instead, the user is advised to use the
1315 * energy group exclusions in the mdp file
1318 "\nThe QM subsystem is divided over multiple topologies. "
1319 "The mutual non-bonded interactions cannot be excluded. "
1320 "There are two ways to achieve this:\n\n"
1321 "1) merge the topologies, such that the atoms of the QM "
1322 "subsystem are all present in one single topology file. "
1323 "In this case this warning will dissappear\n\n"
1324 "2) exclude the non-bonded interactions explicitly via the "
1325 "energygrp-excl option in the mdp file. if this is the case "
1326 "this warning may be ignored"