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->bKeepAllGeneratedDihedrals,r1->bKeepAllGeneratedDihedrals);
893 check_restp_type("nrexcl",r0->nrexcl,r1->nrexcl);
894 check_restp_type("HH14",r0->bGenerateHH14Interactions,r1->bGenerateHH14Interactions);
895 check_restp_type("remove dihedrals",r0->bRemoveDihedralIfWithImproper,r1->bRemoveDihedralIfWithImproper);
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. Fix your terminal residues so that they match the residue database (.rtp) entries, or provide terminal database entries (.tdb).");
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. Fix your coordinate file, add a new terminal database entry (.tdb), or select the proper existing 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,plist,excls,hb,bAllowMissing);
1507 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1508 if (plist[F_CMAP].nr > 0)
1510 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1515 /* set mass of all remaining hydrogen atoms */
1517 do_h_mass(&(plist[F_BONDS]),vsite_type,atoms,mHmult,bDeuterate);
1520 /* Cleanup bonds (sort and rm doubles) */
1521 /* clean_bonds(&(plist[F_BONDS]));*/
1524 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1525 " %4d pairs, %4d bonds and %4d virtual sites\n",
1526 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1527 plist[F_LJ14].nr, plist[F_BONDS].nr,
1528 plist[F_VSITE2].nr +
1529 plist[F_VSITE3].nr +
1530 plist[F_VSITE3FD].nr +
1531 plist[F_VSITE3FAD].nr +
1532 plist[F_VSITE3OUT].nr +
1533 plist[F_VSITE4FD].nr +
1534 plist[F_VSITE4FDN].nr );
1536 print_sums(atoms, FALSE);
1538 if (FALSE == bChargeGroups)
1540 scrub_charge_groups(cgnr, atoms->nr);
1545 for(i=0; i<atoms->nres; i++)
1547 atoms->resinfo[i].nr = i + 1;
1548 atoms->resinfo[i].ic = ' ';
1553 fprintf(stderr,"Writing topology\n");
1554 /* We can copy the bonded types from the first restp,
1555 * since the types have to be identical for all residues in one molecule.
1557 for(i=0; i<ebtsNR; i++) {
1558 bts[i] = restp[0].rb[i].type;
1560 write_top(top_file, posre_fn, molname,
1562 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1566 free_t_hackblock(atoms->nres, &hb);
1567 free_t_restp(atoms->nres, &restp);
1568 gmx_residuetype_destroy(rt);
1570 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1572 for (i=0; i<F_NRE; i++)
1573 sfree(plist[i].param);
1574 for (i=0; i<atoms->nr; i++)