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
51 #include "gmx_fatal.h"
53 #include "gpp_nextnb.h"
66 #include "gen_vsite.h"
69 #include "fflibutil.h"
72 /* this must correspond to enum in pdb2top.h */
73 const char *hh
[ehisNR
] = { "HISD", "HISE", "HISH", "HIS1" };
75 static int missing_atoms(t_restp
*rp
, int resind
,t_atoms
*at
, int i0
, int i
)
79 gmx_bool bFound
, bRet
;
82 for (j
=0; j
<rp
->natom
; j
++)
84 name
=*(rp
->atomname
[j
]);
88 bFound
= (bFound
|| !gmx_strcasecmp(*(at
->atomname
[k
]),name
));
93 fprintf(stderr
,"\nWARNING: "
94 "atom %s is missing in residue %s %d in the pdb file\n",
95 name
,*(at
->resinfo
[resind
].name
),at
->resinfo
[resind
].nr
);
96 if (name
[0]=='H' || name
[0]=='h')
98 fprintf(stderr
," You might need to add atom %s to the hydrogen database of building block %s\n"
99 " in the file %s.hdb (see the manual)\n",
100 name
,*(at
->resinfo
[resind
].rtp
),rp
->filebase
);
102 fprintf(stderr
,"\n");
109 gmx_bool
is_int(double x
)
111 const double tol
= 1e-4;
118 return (fabs(x
-ix
) < tol
);
121 static void swap_strings(char **s
,int i
,int j
)
131 choose_ff(const char *ffsel
,
132 char *forcefield
, int ff_maxlen
,
133 char *ffdir
, int ffdir_maxlen
)
136 char **ffdirs
,**ffs
,**ffs_dir
,*ptr
;
137 int i
,j
,sel
,cwdsel
,nfound
;
138 char buf
[STRLEN
],**desc
;
142 nff
= fflib_search_file_in_dirend(fflib_forcefield_itp(),
143 fflib_forcefield_dir_ext(),
148 gmx_fatal(FARGS
,"No force fields found (files with name '%s' in subdirectories ending on '%s')",
149 fflib_forcefield_itp(),fflib_forcefield_dir_ext());
152 /* Replace with unix path separators */
153 if(DIR_SEPARATOR
!='/')
157 while( (ptr
=strchr(ffdirs
[i
],DIR_SEPARATOR
))!=NULL
)
164 /* Store the force field names in ffs */
169 /* Remove the path from the ffdir name - use our unix standard here! */
170 ptr
= strrchr(ffdirs
[i
],'/');
173 ffs
[i
] = strdup(ffdirs
[i
]);
174 ffs_dir
[i
] = low_gmxlibfn(ffdirs
[i
],FALSE
,FALSE
);
175 if (ffs_dir
[i
] == NULL
)
177 gmx_fatal(FARGS
,"Can no longer find file '%s'",ffdirs
[i
]);
182 ffs
[i
] = strdup(ptr
+1);
183 ffs_dir
[i
] = strdup(ffdirs
[i
]);
185 ffs_dir
[i
][strlen(ffs_dir
[i
])-strlen(ffs
[i
])-1] = '\0';
186 /* Remove the extension from the ffdir name */
187 ffs
[i
][strlen(ffs
[i
])-strlen(fflib_forcefield_dir_ext())] = '\0';
197 if ( strcmp(ffs
[i
],ffsel
)==0 )
199 /* Matching ff name */
203 if( strncmp(ffs_dir
[i
],".",1)==0 )
220 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
221 "current directory. Use interactive selection (not the -ff option) if\n"
222 "you would prefer a different one.\n",ffsel
,nfound
);
227 "Force field '%s' occurs in %d places, but not in the current directory.\n"
228 "Run without the -ff switch and select the force field interactively.",ffsel
,nfound
);
233 gmx_fatal(FARGS
,"Could not find force field '%s' in current directory, install tree or GMXDATA path.",ffsel
);
239 for(i
=0; (i
<nff
); i
++)
241 sprintf(buf
,"%s%c%s%s%c%s",
242 ffs_dir
[i
],DIR_SEPARATOR
,
243 ffs
[i
],fflib_forcefield_dir_ext(),DIR_SEPARATOR
,
244 fflib_forcefield_doc());
247 /* We don't use fflib_open, because we don't want printf's */
248 fp
= ffopen(buf
,"r");
249 snew(desc
[i
],STRLEN
);
250 get_a_line(fp
,desc
[i
],STRLEN
);
255 desc
[i
] = strdup(ffs
[i
]);
258 /* Order force fields from the same dir alphabetically
259 * and put deprecated force fields at the end.
261 for(i
=0; (i
<nff
); i
++)
263 for(j
=i
+1; (j
<nff
); j
++)
265 if (strcmp(ffs_dir
[i
],ffs_dir
[j
]) == 0 &&
266 ((desc
[i
][0] == '[' && desc
[j
][0] != '[') ||
267 ((desc
[i
][0] == '[' || desc
[j
][0] != '[') &&
268 gmx_strcasecmp(desc
[i
],desc
[j
]) > 0)))
270 swap_strings(ffdirs
,i
,j
);
271 swap_strings(ffs
,i
,j
);
272 swap_strings(desc
,i
,j
);
277 printf("\nSelect the Force Field:\n");
278 for(i
=0; (i
<nff
); i
++)
280 if (i
== 0 || strcmp(ffs_dir
[i
-1],ffs_dir
[i
]) != 0)
282 if( strcmp(ffs_dir
[i
],".")==0 )
284 printf("From current directory:\n");
288 printf("From '%s':\n",ffs_dir
[i
]);
291 printf("%2d: %s\n",i
+1,desc
[i
]);
298 pret
= fgets(buf
,STRLEN
,stdin
);
302 sscanf(buf
,"%d",&sel
);
306 while ( pret
==NULL
|| (sel
< 0) || (sel
>= nff
));
308 /* Check for a current limitation of the fflib code.
309 * It will always read from the first ff directory in the list.
310 * This check assumes that the order of ffs matches the order
311 * in which fflib_open searches ff library files.
315 if (strcmp(ffs
[i
],ffs
[sel
]) == 0)
317 gmx_fatal(FARGS
,"Can only select the first of multiple force field entries with directory name '%s%s' in the list. If you want to use the next entry, run pdb2gmx in a different directory or rename or move the force field directory present in the current working directory.",
318 ffs
[sel
],fflib_forcefield_dir_ext());
327 if (strlen(ffs
[sel
]) >= (size_t)ff_maxlen
)
329 gmx_fatal(FARGS
,"Length of force field name (%d) >= maxlen (%d)",
330 strlen(ffs
[sel
]),ff_maxlen
);
332 strcpy(forcefield
,ffs
[sel
]);
334 if (strlen(ffdirs
[sel
]) >= (size_t)ffdir_maxlen
)
336 gmx_fatal(FARGS
,"Length of force field dir (%d) >= maxlen (%d)",
337 strlen(ffdirs
[sel
]),ffdir_maxlen
);
339 strcpy(ffdir
,ffdirs
[sel
]);
341 for(i
=0; (i
<nff
); i
++)
352 void choose_watermodel(const char *wmsel
,const char *ffdir
,
355 const char *fn_watermodels
="watermodels.dat";
356 char fn_list
[STRLEN
];
363 if (strcmp(wmsel
,"none") == 0)
369 else if (strcmp(wmsel
,"select") != 0)
371 *watermodel
= strdup(wmsel
);
376 sprintf(fn_list
,"%s%c%s",ffdir
,DIR_SEPARATOR
,fn_watermodels
);
378 if (!fflib_fexist(fn_list
))
380 fprintf(stderr
,"No file '%s' found, will not include a water model\n",
387 fp
= fflib_open(fn_list
);
388 printf("\nSelect the Water Model:\n");
391 while (get_a_line(fp
,buf
,STRLEN
))
394 snew(model
[nwm
],STRLEN
);
395 sscanf(buf
,"%s%n",model
[nwm
],&i
);
399 fprintf(stderr
,"%2d: %s\n",nwm
+1,buf
+i
);
408 fprintf(stderr
,"%2d: %s\n",nwm
+1,"None");
412 pret
= fgets(buf
,STRLEN
,stdin
);
416 sscanf(buf
,"%d",&sel
);
420 while (pret
== NULL
|| sel
< 0 || sel
> nwm
);
428 *watermodel
= strdup(model
[sel
]);
438 static int name2type(t_atoms
*at
, int **cgnr
, gpp_atomtype_t atype
,
439 t_restp restp
[], gmx_residuetype_t rt
)
441 int i
,j
,prevresind
,resind
,i0
,prevcg
,cg
,curcg
;
443 gmx_bool bProt
, bNterm
;
460 for(i
=0; (i
<at
->nr
); i
++) {
462 if (at
->atom
[i
].resind
!= resind
) {
463 resind
= at
->atom
[i
].resind
;
464 bProt
= gmx_residuetype_is_protein(rt
,*(at
->resinfo
[resind
].name
));
465 bNterm
=bProt
&& (resind
== 0);
467 nmissat
+= missing_atoms(&restp
[prevresind
],prevresind
,at
,i0
,i
);
471 if (at
->atom
[i
].m
== 0) {
473 fprintf(debug
,"atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
474 i
+1,*(at
->atomname
[i
]),curcg
,prevcg
,
475 j
==NOTSET
? NOTSET
: restp
[resind
].cgnr
[j
]);
478 name
=*(at
->atomname
[i
]);
479 j
=search_jtype(&restp
[resind
],name
,bNterm
);
480 at
->atom
[i
].type
= restp
[resind
].atom
[j
].type
;
481 at
->atom
[i
].q
= restp
[resind
].atom
[j
].q
;
482 at
->atom
[i
].m
= get_atomtype_massA(restp
[resind
].atom
[j
].type
,
484 cg
= restp
[resind
].cgnr
[j
];
485 /* A charge group number -1 signals a separate charge group
488 if ( (cg
== -1) || (cg
!= prevcg
) || (resind
!= prevresind
) ) {
493 fprintf(debug
,"atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
494 i
+1,*(at
->atomname
[i
]),curcg
,qt
,is_int(qt
));
503 at
->atom
[i
].typeB
= at
->atom
[i
].type
;
504 at
->atom
[i
].qB
= at
->atom
[i
].q
;
505 at
->atom
[i
].mB
= at
->atom
[i
].m
;
507 nmissat
+= missing_atoms(&restp
[resind
],resind
,at
,i0
,i
);
512 static void print_top_heavy_H(FILE *out
, real mHmult
)
515 fprintf(out
,"; Using deuterium instead of hydrogen\n\n");
516 else if (mHmult
== 4.0)
517 fprintf(out
,"#define HEAVY_H\n\n");
518 else if (mHmult
!= 1.0)
519 fprintf(stderr
,"WARNING: unsupported proton mass multiplier (%g) "
520 "in pdb2top\n",mHmult
);
523 void print_top_comment(FILE *out
,
524 const char *filename
,
525 const char *generator
,
530 char ffdir_parent
[STRLEN
];
533 nice_header(out
,filename
);
534 fprintf(out
,";\tThis is a %s topology file\n;\n",bITP
? "include" : "standalone");
535 fprintf(out
,";\tIt was generated using program:\n;\t%s\n;\n",
536 (NULL
== generator
) ? "unknown" : generator
);
537 fprintf(out
,";\tCommand line was:\n;\t%s\n;\n",command_line());
539 if(strchr(ffdir
,'/')==NULL
)
541 fprintf(out
,";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
543 else if(ffdir
[0]=='.')
545 fprintf(out
,";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
549 strncpy(ffdir_parent
,ffdir
,STRLEN
-1);
550 ffdir_parent
[STRLEN
-1]='\0'; /*make sure it is 0-terminated even for long string*/
551 p
=strrchr(ffdir_parent
,'/');
556 ";\tForce field data was read from:\n"
560 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
561 ";\tforce field must either be present in the current directory, or the location\n"
562 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
567 void print_top_header(FILE *out
,const char *filename
,
568 const char *title
,gmx_bool bITP
,const char *ffdir
,real mHmult
)
572 print_top_comment(out
,filename
,title
,ffdir
,bITP
);
574 print_top_heavy_H(out
, mHmult
);
575 fprintf(out
,"; Include forcefield parameters\n");
577 p
=strrchr(ffdir
,'/');
578 p
= (ffdir
[0]=='.' || p
==NULL
) ? ffdir
: p
+1;
580 fprintf(out
,"#include \"%s/%s\"\n\n",p
,fflib_forcefield_itp());
583 static void print_top_posre(FILE *out
,const char *pr
)
585 fprintf(out
,"; Include Position restraint file\n");
586 fprintf(out
,"#ifdef POSRES\n");
587 fprintf(out
,"#include \"%s\"\n",pr
);
588 fprintf(out
,"#endif\n\n");
591 static void print_top_water(FILE *out
,const char *ffdir
,const char *water
)
596 fprintf(out
,"; Include water topology\n");
598 p
=strrchr(ffdir
,'/');
599 p
= (ffdir
[0]=='.' || p
==NULL
) ? ffdir
: p
+1;
600 fprintf(out
,"#include \"%s/%s.itp\"\n",p
,water
);
603 fprintf(out
,"#ifdef POSRES_WATER\n");
604 fprintf(out
,"; Position restraint for each water oxygen\n");
605 fprintf(out
,"[ position_restraints ]\n");
606 fprintf(out
,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
607 fprintf(out
,"%4d %4d %10g %10g %10g\n",1,1,1000.0,1000.0,1000.0);
608 fprintf(out
,"#endif\n");
611 sprintf(buf
,"%s/ions.itp",p
);
613 if (fflib_fexist(buf
))
615 fprintf(out
,"; Include topology for ions\n");
616 fprintf(out
,"#include \"%s\"\n",buf
);
621 static void print_top_system(FILE *out
, const char *title
)
623 fprintf(out
,"[ %s ]\n",dir2str(d_system
));
624 fprintf(out
,"; Name\n");
625 fprintf(out
,"%s\n\n",title
[0]?title
:"Protein");
628 void print_top_mols(FILE *out
,
629 const char *title
, const char *ffdir
, const char *water
,
630 int nincl
, char **incls
, int nmol
, t_mols
*mols
)
636 fprintf(out
,"; Include chain topologies\n");
637 for (i
=0; (i
<nincl
); i
++) {
638 incl
= strrchr(incls
[i
],DIR_SEPARATOR
);
642 /* Remove the path from the include name */
645 fprintf(out
,"#include \"%s\"\n",incl
);
652 print_top_water(out
,ffdir
,water
);
654 print_top_system(out
, title
);
657 fprintf(out
,"[ %s ]\n",dir2str(d_molecules
));
658 fprintf(out
,"; %-15s %5s\n","Compound","#mols");
659 for (i
=0; (i
<nmol
); i
++)
660 fprintf(out
,"%-15s %5d\n",mols
[i
].name
,mols
[i
].nr
);
664 void write_top(FILE *out
, char *pr
,char *molname
,
665 t_atoms
*at
,gmx_bool bRTPresname
,
666 int bts
[],t_params plist
[],t_excls excls
[],
667 gpp_atomtype_t atype
,int *cgnr
, int nrexcl
)
668 /* NOTE: nrexcl is not the size of *excl! */
670 if (at
&& atype
&& cgnr
) {
671 fprintf(out
,"[ %s ]\n",dir2str(d_moleculetype
));
672 fprintf(out
,"; %-15s %5s\n","Name","nrexcl");
673 fprintf(out
,"%-15s %5d\n\n",molname
?molname
:"Protein",nrexcl
);
675 print_atoms(out
, atype
, at
, cgnr
, bRTPresname
);
676 print_bondeds(out
,at
->nr
,d_bonds
, F_BONDS
, bts
[ebtsBONDS
], plist
);
677 print_bondeds(out
,at
->nr
,d_constraints
,F_CONSTR
, 0, plist
);
678 print_bondeds(out
,at
->nr
,d_constraints
,F_CONSTRNC
, 0, plist
);
679 print_bondeds(out
,at
->nr
,d_pairs
, F_LJ14
, 0, plist
);
680 print_excl(out
,at
->nr
,excls
);
681 print_bondeds(out
,at
->nr
,d_angles
, F_ANGLES
, bts
[ebtsANGLES
],plist
);
682 print_bondeds(out
,at
->nr
,d_dihedrals
, F_PDIHS
, bts
[ebtsPDIHS
], plist
);
683 print_bondeds(out
,at
->nr
,d_dihedrals
, F_IDIHS
, bts
[ebtsIDIHS
], plist
);
684 print_bondeds(out
,at
->nr
,d_cmap
, F_CMAP
, bts
[ebtsCMAP
], plist
);
685 print_bondeds(out
,at
->nr
,d_polarization
,F_POLARIZATION
, 0, plist
);
686 print_bondeds(out
,at
->nr
,d_thole_polarization
,F_THOLE_POL
,0, plist
);
687 print_bondeds(out
,at
->nr
,d_vsites2
, F_VSITE2
, 0, plist
);
688 print_bondeds(out
,at
->nr
,d_vsites3
, F_VSITE3
, 0, plist
);
689 print_bondeds(out
,at
->nr
,d_vsites3
, F_VSITE3FD
, 0, plist
);
690 print_bondeds(out
,at
->nr
,d_vsites3
, F_VSITE3FAD
,0, plist
);
691 print_bondeds(out
,at
->nr
,d_vsites3
, F_VSITE3OUT
,0, plist
);
692 print_bondeds(out
,at
->nr
,d_vsites4
, F_VSITE4FD
, 0, plist
);
693 print_bondeds(out
,at
->nr
,d_vsites4
, F_VSITE4FDN
, 0, plist
);
696 print_top_posre(out
,pr
);
700 static atom_id
search_res_atom(const char *type
,int resind
,
702 const char *bondtype
,gmx_bool bAllowMissing
)
706 for(i
=0; (i
<atoms
->nr
); i
++)
708 if (atoms
->atom
[i
].resind
== resind
)
710 return search_atom(type
,i
,atoms
,bondtype
,bAllowMissing
);
717 static void do_ssbonds(t_params
*ps
,t_atoms
*atoms
,
718 int nssbonds
,t_ssbond
*ssbonds
,gmx_bool bAllowMissing
)
723 for(i
=0; (i
<nssbonds
); i
++) {
724 ri
= ssbonds
[i
].res1
;
725 rj
= ssbonds
[i
].res2
;
726 ai
= search_res_atom(ssbonds
[i
].a1
,ri
,atoms
,
727 "special bond",bAllowMissing
);
728 aj
= search_res_atom(ssbonds
[i
].a2
,rj
,atoms
,
729 "special bond",bAllowMissing
);
730 if ((ai
== NO_ATID
) || (aj
== NO_ATID
))
731 gmx_fatal(FARGS
,"Trying to make impossible special bond (%s-%s)!",
732 ssbonds
[i
].a1
,ssbonds
[i
].a2
);
733 add_param(ps
,ai
,aj
,NULL
,NULL
);
737 static void at2bonds(t_params
*psb
, t_hackblock
*hb
,
740 real long_bond_dist
, real short_bond_dist
,
741 gmx_bool bAllowMissing
)
745 real dist2
, long_bond_dist2
, short_bond_dist2
;
748 long_bond_dist2
= sqr(long_bond_dist
);
749 short_bond_dist2
= sqr(short_bond_dist
);
756 fprintf(stderr
,"Making bonds...\n");
758 for(resind
=0; (resind
< atoms
->nres
) && (i
<atoms
->nr
); resind
++) {
759 /* add bonds from list of bonded interactions */
760 for(j
=0; j
< hb
[resind
].rb
[ebtsBONDS
].nb
; j
++) {
761 /* Unfortunately we can not issue errors or warnings
762 * for missing atoms in bonds, as the hydrogens and terminal atoms
763 * have not been added yet.
765 ai
=search_atom(hb
[resind
].rb
[ebtsBONDS
].b
[j
].AI
,i
,atoms
,
767 aj
=search_atom(hb
[resind
].rb
[ebtsBONDS
].b
[j
].AJ
,i
,atoms
,
769 if (ai
!= NO_ATID
&& aj
!= NO_ATID
) {
770 dist2
= distance2(x
[ai
],x
[aj
]);
771 if (dist2
> long_bond_dist2
)
773 fprintf(stderr
,"Warning: Long Bond (%d-%d = %g nm)\n",
774 ai
+1,aj
+1,sqrt(dist2
));
776 else if (dist2
< short_bond_dist2
)
778 fprintf(stderr
,"Warning: Short Bond (%d-%d = %g nm)\n",
779 ai
+1,aj
+1,sqrt(dist2
));
781 add_param(psb
,ai
,aj
,NULL
,hb
[resind
].rb
[ebtsBONDS
].b
[j
].s
);
784 /* add bonds from list of hacks (each added atom gets a bond) */
785 while( (i
<atoms
->nr
) && (atoms
->atom
[i
].resind
== resind
) ) {
786 for(j
=0; j
< hb
[resind
].nhack
; j
++)
787 if ( ( hb
[resind
].hack
[j
].tp
> 0 ||
788 hb
[resind
].hack
[j
].oname
==NULL
) &&
789 strcmp(hb
[resind
].hack
[j
].AI
,*(atoms
->atomname
[i
])) == 0 ) {
790 switch(hb
[resind
].hack
[j
].tp
) {
791 case 9: /* COOH terminus */
792 add_param(psb
,i
,i
+1,NULL
,NULL
); /* C-O */
793 add_param(psb
,i
,i
+2,NULL
,NULL
); /* C-OA */
794 add_param(psb
,i
+2,i
+3,NULL
,NULL
); /* OA-H */
797 for(k
=0; (k
<hb
[resind
].hack
[j
].nr
); k
++)
798 add_param(psb
,i
,i
+k
+1,NULL
,NULL
);
803 /* we're now at the start of the next residue */
807 static int pcompar(const void *a
, const void *b
)
818 return strlen(pb
->s
) - strlen(pa
->s
);
823 static void clean_bonds(t_params
*ps
)
829 /* swap atomnumbers in bond if first larger than second: */
830 for(i
=0; (i
<ps
->nr
); i
++)
831 if ( ps
->param
[i
].AJ
< ps
->param
[i
].AI
) {
833 ps
->param
[i
].AI
= ps
->param
[i
].AJ
;
838 qsort(ps
->param
,ps
->nr
,(size_t)sizeof(ps
->param
[0]),pcompar
);
840 /* remove doubles, keep the first one always. */
842 for(i
=1; (i
<ps
->nr
); i
++) {
843 if ((ps
->param
[i
].AI
!= ps
->param
[j
-1].AI
) ||
844 (ps
->param
[i
].AJ
!= ps
->param
[j
-1].AJ
) ) {
846 cp_param(&(ps
->param
[j
]),&(ps
->param
[i
]));
851 fprintf(stderr
,"Number of bonds was %d, now %d\n",ps
->nr
,j
);
855 fprintf(stderr
,"No bonds\n");
858 void print_sums(t_atoms
*atoms
, gmx_bool bSystem
)
871 for(i
=0; (i
<atoms
->nr
); i
++) {
873 qtot
+=atoms
->atom
[i
].q
;
876 fprintf(stderr
,"Total mass%s %.3f a.m.u.\n",where
,m
);
877 fprintf(stderr
,"Total charge%s %.3f e\n",where
,qtot
);
880 static void check_restp_type(const char *name
,int t1
,int t2
)
884 gmx_fatal(FARGS
,"Residues in one molecule have a different '%s' type: %d and %d",name
,t1
,t2
);
888 static void check_restp_types(t_restp
*r0
,t_restp
*r1
)
892 check_restp_type("all dihedrals",r0
->bAlldih
,r1
->bAlldih
);
893 check_restp_type("nrexcl",r0
->nrexcl
,r1
->nrexcl
);
894 check_restp_type("HH14",r0
->HH14
,r1
->HH14
);
895 check_restp_type("remove dihedrals",r0
->bRemoveDih
,r1
->bRemoveDih
);
897 for(i
=0; i
<ebtsNR
; i
++)
899 check_restp_type(btsNames
[i
],r0
->rb
[i
].type
,r1
->rb
[i
].type
);
903 void add_atom_to_restp(t_restp
*restp
,int resnr
,int at_start
,const t_hack
*hack
)
907 const char *Hnum
="123456";
911 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
913 *restp->atomname[at_start], resnr, restp->resname);
915 strcpy(buf
, hack
->nname
);
916 buf
[strlen(buf
)+1]='\0';
919 buf
[strlen(buf
)]='-';
922 restp
->natom
+= hack
->nr
;
923 srenew(restp
->atom
, restp
->natom
);
924 srenew(restp
->atomname
, restp
->natom
);
925 srenew(restp
->cgnr
, restp
->natom
);
927 for(k
=restp
->natom
-1; k
> at_start
+hack
->nr
; k
--)
930 restp
->atom
[k
- hack
->nr
];
932 restp
->atomname
[k
- hack
->nr
];
934 restp
->cgnr
[k
- hack
->nr
];
937 for(k
=0; k
< hack
->nr
; k
++)
939 /* set counter in atomname */
942 buf
[strlen(buf
)-1] = Hnum
[k
];
944 snew( restp
->atomname
[at_start
+1+k
], 1);
945 restp
->atom
[at_start
+1+k
] = *hack
->atom
;
946 *restp
->atomname
[at_start
+1+k
] = strdup(buf
);
947 if ( hack
->cgnr
!= NOTSET
)
949 restp
->cgnr
[at_start
+1+k
] = hack
->cgnr
;
953 restp
->cgnr
[at_start
+1+k
] = restp
->cgnr
[at_start
];
958 void get_hackblocks_rtp(t_hackblock
**hb
, t_restp
**restp
,
959 int nrtp
, t_restp rtp
[],
960 int nres
, t_resinfo
*resinfo
,
962 t_hackblock
**ntdb
, t_hackblock
**ctdb
,
969 const char *Hnum
="123456";
975 /* first the termini */
976 for(i
=0; i
<nterpairs
; i
++) {
977 if (rn
[i
] >= 0 && ntdb
[i
] != NULL
) {
978 copy_t_hackblock(ntdb
[i
], &(*hb
)[rn
[i
]]);
980 if (rc
[i
] >= 0 && ctdb
[i
] != NULL
) {
981 merge_t_hackblock(ctdb
[i
], &(*hb
)[rc
[i
]]);
985 /* then the whole rtp */
986 for(i
=0; i
< nres
; i
++) {
987 /* Here we allow a mismatch of one character when looking for the rtp entry.
988 * For such a mismatch there should be only one mismatching name.
989 * This is mainly useful for small molecules such as ions.
990 * Note that this will usually not work for protein, DNA and RNA,
991 * since there the residue names should be listed in residuetypes.dat
992 * and an error will have been generated earlier in the process.
994 key
= *resinfo
[i
].rtp
;
995 snew(resinfo
[i
].rtp
,1);
996 *resinfo
[i
].rtp
= search_rtp(key
,nrtp
,rtp
);
997 res
= get_restp(*resinfo
[i
].rtp
,nrtp
,rtp
);
998 copy_t_restp(res
, &(*restp
)[i
]);
1000 /* Check that we do not have different bonded types in one molecule */
1001 check_restp_types(&(*restp
)[0],&(*restp
)[i
]);
1004 for(j
=0; j
<nterpairs
&& tern
==-1; j
++) {
1010 for(j
=0; j
<nterpairs
&& terc
== -1; j
++) {
1015 bRM
= merge_t_bondeds(res
->rb
, (*hb
)[i
].rb
,tern
>=0,terc
>=0);
1017 if (bRM
&& ((tern
>= 0 && ntdb
[tern
] == NULL
) ||
1018 (terc
>= 0 && ctdb
[terc
] == NULL
))) {
1019 gmx_fatal(FARGS
,"There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Edit a .n.tdb and/or .c.tdb file.");
1021 if (bRM
&& ((tern
>= 0 && ntdb
[tern
]->nhack
== 0) ||
1022 (terc
>= 0 && ctdb
[terc
]->nhack
== 0))) {
1023 gmx_fatal(FARGS
,"There is a dangling bond at at least one of the terminal ends. Select a proper terminal entry.");
1027 /* now perform t_hack's on t_restp's,
1028 i.e. add's and deletes from termini database will be
1029 added to/removed from residue topology
1030 we'll do this on one big dirty loop, so it won't make easy reading! */
1031 for(i
=0; i
< nres
; i
++)
1033 for(j
=0; j
< (*hb
)[i
].nhack
; j
++)
1035 if ( (*hb
)[i
].hack
[j
].nr
)
1037 /* find atom in restp */
1038 for(l
=0; l
< (*restp
)[i
].natom
; l
++)
1039 if ( ( (*hb
)[i
].hack
[j
].oname
==NULL
&&
1040 strcmp((*hb
)[i
].hack
[j
].AI
, *(*restp
)[i
].atomname
[l
])==0 ) ||
1041 ( (*hb
)[i
].hack
[j
].oname
!=NULL
&&
1042 strcmp((*hb
)[i
].hack
[j
].oname
,*(*restp
)[i
].atomname
[l
])==0 ) )
1044 if (l
== (*restp
)[i
].natom
)
1046 /* If we are doing an atom rename only, we don't need
1047 * to generate a fatal error if the old name is not found
1050 /* Deleting can happen also only on the input atoms,
1051 * not necessarily always on the rtp entry.
1053 if (!((*hb
)[i
].hack
[j
].oname
!= NULL
&&
1054 (*hb
)[i
].hack
[j
].nname
!= NULL
) &&
1055 !((*hb
)[i
].hack
[j
].oname
!= NULL
&&
1056 (*hb
)[i
].hack
[j
].nname
== NULL
))
1059 "atom %s not found in buiding block %d%s "
1060 "while combining tdb and rtp",
1061 (*hb
)[i
].hack
[j
].oname
!=NULL
?
1062 (*hb
)[i
].hack
[j
].oname
: (*hb
)[i
].hack
[j
].AI
,
1063 i
+1,*resinfo
[i
].rtp
);
1068 if ( (*hb
)[i
].hack
[j
].oname
==NULL
) {
1070 add_atom_to_restp(&(*restp
)[i
],resinfo
[i
].nr
,l
,
1076 if ( (*hb
)[i
].hack
[j
].nname
==NULL
) {
1077 /* we're deleting */
1079 fprintf(debug
, "deleting atom %s from res %d%s in rtp\n",
1080 *(*restp
)[i
].atomname
[l
],
1081 i
+1,(*restp
)[i
].resname
);
1082 /* shift the rest */
1083 (*restp
)[i
].natom
--;
1084 for(k
=l
; k
< (*restp
)[i
].natom
; k
++) {
1085 (*restp
)[i
].atom
[k
] = (*restp
)[i
].atom
[k
+1];
1086 (*restp
)[i
].atomname
[k
] = (*restp
)[i
].atomname
[k
+1];
1087 (*restp
)[i
].cgnr
[k
] = (*restp
)[i
].cgnr
[k
+1];
1089 /* give back space */
1090 srenew((*restp
)[i
].atom
, (*restp
)[i
].natom
);
1091 srenew((*restp
)[i
].atomname
, (*restp
)[i
].natom
);
1092 srenew((*restp
)[i
].cgnr
, (*restp
)[i
].natom
);
1093 } else { /* nname != NULL */
1094 /* we're replacing */
1096 fprintf(debug
, "replacing atom %s by %s in res %d%s in rtp\n",
1097 *(*restp
)[i
].atomname
[l
], (*hb
)[i
].hack
[j
].nname
,
1098 i
+1,(*restp
)[i
].resname
);
1099 snew( (*restp
)[i
].atomname
[l
], 1);
1100 (*restp
)[i
].atom
[l
] = *(*hb
)[i
].hack
[j
].atom
;
1101 *(*restp
)[i
].atomname
[l
] = strdup((*hb
)[i
].hack
[j
].nname
);
1102 if ( (*hb
)[i
].hack
[j
].cgnr
!= NOTSET
)
1103 (*restp
)[i
].cgnr
[l
] = (*hb
)[i
].hack
[j
].cgnr
;
1112 static gmx_bool
atomname_cmp_nr(const char *anm
,t_hack
*hack
,int *nr
)
1119 return (gmx_strcasecmp(anm
,hack
->nname
) == 0);
1123 if (isdigit(anm
[strlen(anm
)-1]))
1125 *nr
= anm
[strlen(anm
)-1] - '0';
1131 if (*nr
<= 0 || *nr
> hack
->nr
)
1137 return (strlen(anm
) == strlen(hack
->nname
) + 1 &&
1138 gmx_strncasecmp(anm
,hack
->nname
,strlen(hack
->nname
)) == 0);
1143 static gmx_bool
match_atomnames_with_rtp_atom(t_atoms
*pdba
,rvec
*x
,int atind
,
1144 t_restp
*rptr
,t_hackblock
*hbr
,
1151 char *start_at
,buf
[STRLEN
];
1153 gmx_bool bReplaceReplace
,bFoundInAdd
;
1156 oldnm
= *pdba
->atomname
[atind
];
1157 resnr
= pdba
->resinfo
[pdba
->atom
[atind
].resind
].nr
;
1160 for(j
=0; j
<hbr
->nhack
; j
++)
1162 if (hbr
->hack
[j
].oname
!= NULL
&& hbr
->hack
[j
].nname
!= NULL
&&
1163 gmx_strcasecmp(oldnm
,hbr
->hack
[j
].oname
) == 0)
1165 /* This is a replace entry. */
1166 /* Check if we are not replacing a replaced atom. */
1167 bReplaceReplace
= FALSE
;
1168 for(k
=0; k
<hbr
->nhack
; k
++) {
1170 hbr
->hack
[k
].oname
!= NULL
&& hbr
->hack
[k
].nname
!= NULL
&&
1171 gmx_strcasecmp(hbr
->hack
[k
].nname
,hbr
->hack
[j
].oname
) == 0)
1173 /* The replace in hack[j] replaces an atom that
1174 * was already replaced in hack[k], we do not want
1175 * second or higher level replaces at this stage.
1177 bReplaceReplace
= TRUE
;
1180 if (bReplaceReplace
)
1182 /* Skip this replace. */
1186 /* This atom still has the old name, rename it */
1187 newnm
= hbr
->hack
[j
].nname
;
1188 for(k
=0; k
<rptr
->natom
; k
++)
1190 if (gmx_strcasecmp(newnm
,*rptr
->atomname
[k
]) == 0)
1195 if (k
== rptr
->natom
)
1197 /* The new name is not present in the rtp.
1198 * We need to apply the replace also to the rtp entry.
1201 /* We need to find the add hack that can add this atom
1202 * to find out after which atom it should be added.
1204 bFoundInAdd
= FALSE
;
1205 for(k
=0; k
<hbr
->nhack
; k
++)
1207 if (hbr
->hack
[k
].oname
== NULL
&&
1208 hbr
->hack
[k
].nname
!= NULL
&&
1209 atomname_cmp_nr(newnm
,&hbr
->hack
[k
],&anmnr
))
1213 start_at
= hbr
->hack
[k
].a
[0];
1217 sprintf(buf
,"%s%d",hbr
->hack
[k
].nname
,anmnr
-1);
1220 for(start_nr
=0; start_nr
<rptr
->natom
; start_nr
++)
1222 if (gmx_strcasecmp(start_at
,(*rptr
->atomname
[start_nr
])) == 0)
1227 if (start_nr
== rptr
->natom
)
1229 gmx_fatal(FARGS
,"Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1230 start_at
,rptr
->resname
,newnm
);
1232 /* We can add the atom after atom start_nr */
1233 add_atom_to_restp(rptr
,resnr
,start_nr
,
1242 gmx_fatal(FARGS
,"Could not find an 'add' entry for atom named '%s' corresponding to the 'replace' entry from atom name '%s' to '%s' for tdb or hdb database of residue type '%s'",
1244 hbr
->hack
[j
].oname
,hbr
->hack
[j
].nname
,
1251 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1252 oldnm
,rptr
->resname
,resnr
,newnm
);
1254 /* Rename the atom in pdba */
1255 snew(pdba
->atomname
[atind
],1);
1256 *pdba
->atomname
[atind
] = strdup(newnm
);
1258 else if (hbr
->hack
[j
].oname
!= NULL
&& hbr
->hack
[j
].nname
== NULL
&&
1259 gmx_strcasecmp(oldnm
,hbr
->hack
[j
].oname
) == 0)
1261 /* This is a delete entry, check if this atom is present
1262 * in the rtp entry of this residue.
1264 for(k
=0; k
<rptr
->natom
; k
++)
1266 if (gmx_strcasecmp(oldnm
,*rptr
->atomname
[k
]) == 0)
1271 if (k
== rptr
->natom
)
1273 /* This atom is not present in the rtp entry,
1274 * delete is now from the input pdba.
1278 printf("Deleting atom '%s' in residue '%s' %d\n",
1279 oldnm
,rptr
->resname
,resnr
);
1281 /* We should free the atom name,
1282 * but it might be used multiple times in the symtab.
1283 * sfree(pdba->atomname[atind]);
1285 for(k
=atind
+1; k
<pdba
->nr
; k
++)
1287 pdba
->atom
[k
-1] = pdba
->atom
[k
];
1288 pdba
->atomname
[k
-1] = pdba
->atomname
[k
];
1289 copy_rvec(x
[k
],x
[k
-1]);
1300 void match_atomnames_with_rtp(t_restp restp
[],t_hackblock hb
[],
1301 t_atoms
*pdba
,rvec
*x
,
1310 char *start_at
,buf
[STRLEN
];
1312 gmx_bool bFoundInAdd
;
1314 for(i
=0; i
<pdba
->nr
; i
++)
1316 oldnm
= *pdba
->atomname
[i
];
1317 resnr
= pdba
->resinfo
[pdba
->atom
[i
].resind
].nr
;
1318 rptr
= &restp
[pdba
->atom
[i
].resind
];
1319 for(j
=0; (j
<rptr
->natom
); j
++)
1321 if (gmx_strcasecmp(oldnm
,*(rptr
->atomname
[j
])) == 0)
1326 if (j
== rptr
->natom
)
1328 /* Not found yet, check if we have to rename this atom */
1329 if (match_atomnames_with_rtp_atom(pdba
,x
,i
,
1330 rptr
,&(hb
[pdba
->atom
[i
].resind
]),
1333 /* We deleted this atom, decrease the atom counter by 1. */
1340 #define NUM_CMAP_ATOMS 5
1341 static void gen_cmap(t_params
*psb
, t_restp
*restp
, t_atoms
*atoms
, gmx_residuetype_t rt
)
1345 t_resinfo
*resinfo
= atoms
->resinfo
;
1346 int nres
= atoms
->nres
;
1348 atom_id cmap_atomid
[NUM_CMAP_ATOMS
];
1349 int cmap_chainnum
=-1, this_residue_index
;
1356 fprintf(stderr
,"Making cmap torsions...");
1358 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1359 * therefore we get a valgrind invalid 4 byte read error with atom am */
1360 for(residx
=0; residx
<nres
-1; residx
++)
1362 /* Add CMAP terms from the list of CMAP interactions */
1363 for(j
=0;j
<restp
[residx
].rb
[ebtsCMAP
].nb
; j
++)
1366 /* Loop over atoms in a candidate CMAP interaction and
1367 * check that they exist, are from the same chain and are
1368 * from residues labelled as protein. */
1369 for(k
= 0; k
< NUM_CMAP_ATOMS
&& bAddCMAP
; k
++)
1371 cmap_atomid
[k
] = search_atom(restp
[residx
].rb
[ebtsCMAP
].b
[j
].a
[k
],
1373 bAddCMAP
= bAddCMAP
&& (cmap_atomid
[k
] != NO_ATID
);
1376 /* This break is necessary, because cmap_atomid[k]
1377 * == NO_ATID cannot be safely used as an index
1378 * into the atom array. */
1381 this_residue_index
= atoms
->atom
[cmap_atomid
[k
]].resind
;
1384 cmap_chainnum
= resinfo
[this_residue_index
].chainnum
;
1388 /* Does the residue for this atom have the same
1389 * chain number as the residues for previous
1391 bAddCMAP
= bAddCMAP
&&
1392 cmap_chainnum
== resinfo
[this_residue_index
].chainnum
;
1394 bAddCMAP
= bAddCMAP
&& gmx_residuetype_is_protein(rt
,*(resinfo
[this_residue_index
].name
));
1399 add_cmap_param(psb
,cmap_atomid
[0],cmap_atomid
[1],cmap_atomid
[2],cmap_atomid
[3],cmap_atomid
[4],restp
[residx
].rb
[ebtsCMAP
].b
[j
].s
);
1405 while(atoms
->atom
[i
].resind
<residx
+1)
1412 /* Start the next residue */
1416 scrub_charge_groups(int *cgnr
, int natoms
)
1420 for(i
=0;i
<natoms
;i
++)
1427 void pdb2top(FILE *top_file
, char *posre_fn
, char *molname
,
1428 t_atoms
*atoms
, rvec
**x
, gpp_atomtype_t atype
, t_symtab
*tab
,
1429 int nrtp
, t_restp rtp
[],
1430 t_restp
*restp
, t_hackblock
*hb
,
1431 int nterpairs
,t_hackblock
**ntdb
, t_hackblock
**ctdb
,
1432 gmx_bool bAllowMissing
,
1433 gmx_bool bVsites
, gmx_bool bVsiteAromatics
,
1434 const char *ff
, const char *ffdir
,
1436 int nssbonds
, t_ssbond
*ssbonds
,
1437 real long_bond_dist
, real short_bond_dist
,
1438 gmx_bool bDeuterate
, gmx_bool bChargeGroups
, gmx_bool bCmap
,
1439 gmx_bool bRenumRes
,gmx_bool bRTPresname
)
1445 t_params plist
[F_NRE
];
1452 gmx_residuetype_t rt
;
1455 gmx_residuetype_init(&rt
);
1458 print_resall(debug
, atoms
->nres
, restp
, atype
);
1459 dump_hb(debug
, atoms
->nres
, hb
);
1463 at2bonds(&(plist
[F_BONDS
]), hb
,
1465 long_bond_dist
, short_bond_dist
, bAllowMissing
);
1467 /* specbonds: disulphide bonds & heme-his */
1468 do_ssbonds(&(plist
[F_BONDS
]),
1469 atoms
, nssbonds
, ssbonds
,
1472 nmissat
= name2type(atoms
, &cgnr
, atype
, restp
, rt
);
1475 fprintf(stderr
,"There were %d missing atoms in molecule %s\n",
1478 gmx_fatal(FARGS
,"There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1482 /* Cleanup bonds (sort and rm doubles) */
1483 clean_bonds(&(plist
[F_BONDS
]));
1485 snew(vsite_type
,atoms
->nr
);
1486 for(i
=0; i
<atoms
->nr
; i
++)
1487 vsite_type
[i
]=NOTSET
;
1489 /* determine which atoms will be vsites and add dummy masses
1490 also renumber atom numbers in plist[0..F_NRE]! */
1491 do_vsites(nrtp
, rtp
, atype
, atoms
, tab
, x
, plist
,
1492 &vsite_type
, &cgnr
, mHmult
, bVsiteAromatics
, ffdir
);
1495 /* Make Angles and Dihedrals */
1496 fprintf(stderr
,"Generating angles, dihedrals and pairs...\n");
1497 snew(excls
,atoms
->nr
);
1498 init_nnb(&nnb
,atoms
->nr
,4);
1499 gen_nnb(&nnb
,plist
);
1500 print_nnb(&nnb
,"NNB");
1501 gen_pad(&nnb
,atoms
,restp
[0].nrexcl
,restp
[0].HH14
,
1502 plist
,excls
,hb
,restp
[0].bAlldih
,restp
[0].bRemoveDih
,
1509 gen_cmap(&(plist
[F_CMAP
]), restp
, atoms
, rt
);
1510 if (plist
[F_CMAP
].nr
> 0)
1512 fprintf(stderr
, "There are %4d cmap torsion pairs\n",
1517 /* set mass of all remaining hydrogen atoms */
1519 do_h_mass(&(plist
[F_BONDS
]),vsite_type
,atoms
,mHmult
,bDeuterate
);
1522 /* Cleanup bonds (sort and rm doubles) */
1523 /* clean_bonds(&(plist[F_BONDS]));*/
1526 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1527 " %4d pairs, %4d bonds and %4d virtual sites\n",
1528 plist
[F_PDIHS
].nr
, plist
[F_IDIHS
].nr
, plist
[F_ANGLES
].nr
,
1529 plist
[F_LJ14
].nr
, plist
[F_BONDS
].nr
,
1530 plist
[F_VSITE2
].nr
+
1531 plist
[F_VSITE3
].nr
+
1532 plist
[F_VSITE3FD
].nr
+
1533 plist
[F_VSITE3FAD
].nr
+
1534 plist
[F_VSITE3OUT
].nr
+
1535 plist
[F_VSITE4FD
].nr
+
1536 plist
[F_VSITE4FDN
].nr
);
1538 print_sums(atoms
, FALSE
);
1540 if (FALSE
== bChargeGroups
)
1542 scrub_charge_groups(cgnr
, atoms
->nr
);
1547 for(i
=0; i
<atoms
->nres
; i
++)
1549 atoms
->resinfo
[i
].nr
= i
+ 1;
1550 atoms
->resinfo
[i
].ic
= ' ';
1555 fprintf(stderr
,"Writing topology\n");
1556 /* We can copy the bonded types from the first restp,
1557 * since the types have to be identical for all residues in one molecule.
1559 for(i
=0; i
<ebtsNR
; i
++) {
1560 bts
[i
] = restp
[0].rb
[i
].type
;
1562 write_top(top_file
, posre_fn
, molname
,
1564 bts
, plist
, excls
, atype
, cgnr
, restp
[0].nrexcl
);
1568 free_t_hackblock(atoms
->nres
, &hb
);
1569 free_t_restp(atoms
->nres
, &restp
);
1570 gmx_residuetype_destroy(rt
);
1572 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1574 for (i
=0; i
<F_NRE
; i
++)
1575 sfree(plist
[i
].param
);
1576 for (i
=0; i
<atoms
->nr
; i
++)