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
44 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
57 #include "gmx_fatal.h"
59 #include "gpp_nextnb.h"
72 #include "gen_vsite.h"
75 #include "fflibutil.h"
78 /* this must correspond to enum in pdb2top.h */
79 const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
81 static int missing_atoms(t_restp *rp, int resind,t_atoms *at, int i0, int i)
85 gmx_bool bFound, bRet;
88 for (j=0; j<rp->natom; j++)
90 name=*(rp->atomname[j]);
94 bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]),name));
99 fprintf(stderr,"\nWARNING: "
100 "atom %s is missing in residue %s %d in the pdb file\n",
101 name,*(at->resinfo[resind].name),at->resinfo[resind].nr);
102 if (name[0]=='H' || name[0]=='h')
104 fprintf(stderr," You might need to add atom %s to the hydrogen database of building block %s\n"
105 " in the file %s.hdb (see the manual)\n",
106 name,*(at->resinfo[resind].rtp),rp->filebase);
108 fprintf(stderr,"\n");
115 gmx_bool is_int(double x)
117 const double tol = 1e-4;
124 return (fabs(x-ix) < tol);
127 static void swap_strings(char **s,int i,int j)
137 choose_ff(const char *ffsel,
138 char *forcefield, int ff_maxlen,
139 char *ffdir, int ffdir_maxlen)
142 char **ffdirs,**ffs,**ffs_dir,*ptr;
143 int i,j,sel,cwdsel,nfound;
144 char buf[STRLEN],**desc;
148 nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
149 fflib_forcefield_dir_ext(),
154 gmx_fatal(FARGS,"No force fields found (files with name '%s' in subdirectories ending on '%s')",
155 fflib_forcefield_itp(),fflib_forcefield_dir_ext());
158 /* Replace with unix path separators */
159 if(DIR_SEPARATOR!='/')
163 while( (ptr=strchr(ffdirs[i],DIR_SEPARATOR))!=NULL )
170 /* Store the force field names in ffs */
175 /* Remove the path from the ffdir name - use our unix standard here! */
176 ptr = strrchr(ffdirs[i],'/');
179 ffs[i] = strdup(ffdirs[i]);
180 ffs_dir[i] = low_gmxlibfn(ffdirs[i],FALSE,FALSE);
181 if (ffs_dir[i] == NULL)
183 gmx_fatal(FARGS,"Can no longer find file '%s'",ffdirs[i]);
188 ffs[i] = strdup(ptr+1);
189 ffs_dir[i] = strdup(ffdirs[i]);
191 ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
192 /* Remove the extension from the ffdir name */
193 ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
203 if ( strcmp(ffs[i],ffsel)==0 )
205 /* Matching ff name */
209 if( strncmp(ffs_dir[i],".",1)==0 )
226 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
227 "current directory. Use interactive selection (not the -ff option) if\n"
228 "you would prefer a different one.\n",ffsel,nfound);
233 "Force field '%s' occurs in %d places, but not in the current directory.\n"
234 "Run without the -ff switch and select the force field interactively.",ffsel,nfound);
239 gmx_fatal(FARGS,"Could not find force field '%s' in current directory, install tree or GMXDATA path.",ffsel);
245 for(i=0; (i<nff); i++)
247 sprintf(buf,"%s%c%s%s%c%s",
248 ffs_dir[i],DIR_SEPARATOR,
249 ffs[i],fflib_forcefield_dir_ext(),DIR_SEPARATOR,
250 fflib_forcefield_doc());
253 /* We don't use fflib_open, because we don't want printf's */
254 fp = ffopen(buf,"r");
255 snew(desc[i],STRLEN);
256 get_a_line(fp,desc[i],STRLEN);
261 desc[i] = strdup(ffs[i]);
264 /* Order force fields from the same dir alphabetically
265 * and put deprecated force fields at the end.
267 for(i=0; (i<nff); i++)
269 for(j=i+1; (j<nff); j++)
271 if (strcmp(ffs_dir[i],ffs_dir[j]) == 0 &&
272 ((desc[i][0] == '[' && desc[j][0] != '[') ||
273 ((desc[i][0] == '[' || desc[j][0] != '[') &&
274 gmx_strcasecmp(desc[i],desc[j]) > 0)))
276 swap_strings(ffdirs,i,j);
277 swap_strings(ffs ,i,j);
278 swap_strings(desc ,i,j);
283 printf("\nSelect the Force Field:\n");
284 for(i=0; (i<nff); i++)
286 if (i == 0 || strcmp(ffs_dir[i-1],ffs_dir[i]) != 0)
288 if( strcmp(ffs_dir[i],".")==0 )
290 printf("From current directory:\n");
294 printf("From '%s':\n",ffs_dir[i]);
297 printf("%2d: %s\n",i+1,desc[i]);
304 pret = fgets(buf,STRLEN,stdin);
308 sscanf(buf,"%d",&sel);
312 while ( pret==NULL || (sel < 0) || (sel >= nff));
314 /* Check for a current limitation of the fflib code.
315 * It will always read from the first ff directory in the list.
316 * This check assumes that the order of ffs matches the order
317 * in which fflib_open searches ff library files.
321 if (strcmp(ffs[i],ffs[sel]) == 0)
323 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.",
324 ffs[sel],fflib_forcefield_dir_ext());
333 if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
335 gmx_fatal(FARGS,"Length of force field name (%d) >= maxlen (%d)",
336 strlen(ffs[sel]),ff_maxlen);
338 strcpy(forcefield,ffs[sel]);
340 if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
342 gmx_fatal(FARGS,"Length of force field dir (%d) >= maxlen (%d)",
343 strlen(ffdirs[sel]),ffdir_maxlen);
345 strcpy(ffdir,ffdirs[sel]);
347 for(i=0; (i<nff); i++)
358 void choose_watermodel(const char *wmsel,const char *ffdir,
361 const char *fn_watermodels="watermodels.dat";
362 char fn_list[STRLEN];
369 if (strcmp(wmsel,"none") == 0)
375 else if (strcmp(wmsel,"select") != 0)
377 *watermodel = strdup(wmsel);
382 sprintf(fn_list,"%s%c%s",ffdir,DIR_SEPARATOR,fn_watermodels);
384 if (!fflib_fexist(fn_list))
386 fprintf(stderr,"No file '%s' found, will not include a water model\n",
393 fp = fflib_open(fn_list);
394 printf("\nSelect the Water Model:\n");
397 while (get_a_line(fp,buf,STRLEN))
400 snew(model[nwm],STRLEN);
401 sscanf(buf,"%s%n",model[nwm],&i);
405 fprintf(stderr,"%2d: %s\n",nwm+1,buf+i);
414 fprintf(stderr,"%2d: %s\n",nwm+1,"None");
418 pret = fgets(buf,STRLEN,stdin);
422 sscanf(buf,"%d",&sel);
426 while (pret == NULL || sel < 0 || sel > nwm);
434 *watermodel = strdup(model[sel]);
444 static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
445 t_restp restp[], gmx_residuetype_t rt)
447 int i,j,prevresind,resind,i0,prevcg,cg,curcg;
449 gmx_bool bProt, bNterm;
466 for(i=0; (i<at->nr); i++) {
468 if (at->atom[i].resind != resind) {
469 resind = at->atom[i].resind;
470 bProt = gmx_residuetype_is_protein(rt,*(at->resinfo[resind].name));
471 bNterm=bProt && (resind == 0);
473 nmissat += missing_atoms(&restp[prevresind],prevresind,at,i0,i);
477 if (at->atom[i].m == 0) {
479 fprintf(debug,"atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
480 i+1,*(at->atomname[i]),curcg,prevcg,
481 j==NOTSET ? NOTSET : restp[resind].cgnr[j]);
484 name=*(at->atomname[i]);
485 j=search_jtype(&restp[resind],name,bNterm);
486 at->atom[i].type = restp[resind].atom[j].type;
487 at->atom[i].q = restp[resind].atom[j].q;
488 at->atom[i].m = get_atomtype_massA(restp[resind].atom[j].type,
490 cg = restp[resind].cgnr[j];
491 /* A charge group number -1 signals a separate charge group
494 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) ) {
499 fprintf(debug,"atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
500 i+1,*(at->atomname[i]),curcg,qt,is_int(qt));
509 at->atom[i].typeB = at->atom[i].type;
510 at->atom[i].qB = at->atom[i].q;
511 at->atom[i].mB = at->atom[i].m;
513 nmissat += missing_atoms(&restp[resind],resind,at,i0,i);
518 static void print_top_heavy_H(FILE *out, real mHmult)
521 fprintf(out,"; Using deuterium instead of hydrogen\n\n");
522 else if (mHmult == 4.0)
523 fprintf(out,"#define HEAVY_H\n\n");
524 else if (mHmult != 1.0)
525 fprintf(stderr,"WARNING: unsupported proton mass multiplier (%g) "
526 "in pdb2top\n",mHmult);
529 void print_top_comment(FILE *out,
530 const char *filename,
531 const char *generator,
536 char ffdir_parent[STRLEN];
539 nice_header(out,filename);
540 fprintf(out,";\tThis is a %s topology file\n;\n",bITP ? "include" : "standalone");
541 fprintf(out,";\tIt was generated using program:\n;\t%s\n;\n",
542 (NULL == generator) ? "unknown" : generator);
543 fprintf(out,";\tCommand line was:\n;\t%s\n;\n",command_line());
545 if(strchr(ffdir,'/')==NULL)
547 fprintf(out,";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
549 else if(ffdir[0]=='.')
551 fprintf(out,";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
555 strncpy(ffdir_parent,ffdir,STRLEN-1);
556 ffdir_parent[STRLEN-1]='\0'; /*make sure it is 0-terminated even for long string*/
557 p=strrchr(ffdir_parent,'/');
562 ";\tForce field data was read from:\n"
566 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
567 ";\tforce field must either be present in the current directory, or the location\n"
568 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
573 void print_top_header(FILE *out,const char *filename,
574 const char *title,gmx_bool bITP,const char *ffdir,real mHmult)
578 print_top_comment(out,filename,title,ffdir,bITP);
580 print_top_heavy_H(out, mHmult);
581 fprintf(out,"; Include forcefield parameters\n");
583 p=strrchr(ffdir,'/');
584 p = (ffdir[0]=='.' || p==NULL) ? ffdir : p+1;
586 fprintf(out,"#include \"%s/%s\"\n\n",p,fflib_forcefield_itp());
589 static void print_top_posre(FILE *out,const char *pr)
591 fprintf(out,"; Include Position restraint file\n");
592 fprintf(out,"#ifdef POSRES\n");
593 fprintf(out,"#include \"%s\"\n",pr);
594 fprintf(out,"#endif\n\n");
597 static void print_top_water(FILE *out,const char *ffdir,const char *water)
602 fprintf(out,"; Include water topology\n");
604 p=strrchr(ffdir,'/');
605 p = (ffdir[0]=='.' || p==NULL) ? ffdir : p+1;
606 fprintf(out,"#include \"%s/%s.itp\"\n",p,water);
609 fprintf(out,"#ifdef POSRES_WATER\n");
610 fprintf(out,"; Position restraint for each water oxygen\n");
611 fprintf(out,"[ position_restraints ]\n");
612 fprintf(out,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
613 fprintf(out,"%4d %4d %10g %10g %10g\n",1,1,1000.0,1000.0,1000.0);
614 fprintf(out,"#endif\n");
617 sprintf(buf,"%s/ions.itp",p);
619 if (fflib_fexist(buf))
621 fprintf(out,"; Include topology for ions\n");
622 fprintf(out,"#include \"%s\"\n",buf);
627 static void print_top_system(FILE *out, const char *title)
629 fprintf(out,"[ %s ]\n",dir2str(d_system));
630 fprintf(out,"; Name\n");
631 fprintf(out,"%s\n\n",title[0]?title:"Protein");
634 void print_top_mols(FILE *out,
635 const char *title, const char *ffdir, const char *water,
636 int nincl, char **incls, int nmol, t_mols *mols)
642 fprintf(out,"; Include chain topologies\n");
643 for (i=0; (i<nincl); i++) {
644 incl = strrchr(incls[i],DIR_SEPARATOR);
648 /* Remove the path from the include name */
651 fprintf(out,"#include \"%s\"\n",incl);
658 print_top_water(out,ffdir,water);
660 print_top_system(out, title);
663 fprintf(out,"[ %s ]\n",dir2str(d_molecules));
664 fprintf(out,"; %-15s %5s\n","Compound","#mols");
665 for (i=0; (i<nmol); i++)
666 fprintf(out,"%-15s %5d\n",mols[i].name,mols[i].nr);
670 void write_top(FILE *out, char *pr,char *molname,
671 t_atoms *at,gmx_bool bRTPresname,
672 int bts[],t_params plist[],t_excls excls[],
673 gpp_atomtype_t atype,int *cgnr, int nrexcl)
674 /* NOTE: nrexcl is not the size of *excl! */
676 if (at && atype && cgnr) {
677 fprintf(out,"[ %s ]\n",dir2str(d_moleculetype));
678 fprintf(out,"; %-15s %5s\n","Name","nrexcl");
679 fprintf(out,"%-15s %5d\n\n",molname?molname:"Protein",nrexcl);
681 print_atoms(out, atype, at, cgnr, bRTPresname);
682 print_bondeds(out,at->nr,d_bonds, F_BONDS, bts[ebtsBONDS], plist);
683 print_bondeds(out,at->nr,d_constraints,F_CONSTR, 0, plist);
684 print_bondeds(out,at->nr,d_constraints,F_CONSTRNC, 0, plist);
685 print_bondeds(out,at->nr,d_pairs, F_LJ14, 0, plist);
686 print_excl(out,at->nr,excls);
687 print_bondeds(out,at->nr,d_angles, F_ANGLES, bts[ebtsANGLES],plist);
688 print_bondeds(out,at->nr,d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
689 print_bondeds(out,at->nr,d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
690 print_bondeds(out,at->nr,d_cmap, F_CMAP, bts[ebtsCMAP], plist);
691 print_bondeds(out,at->nr,d_polarization,F_POLARIZATION, 0, plist);
692 print_bondeds(out,at->nr,d_thole_polarization,F_THOLE_POL,0, plist);
693 print_bondeds(out,at->nr,d_vsites2, F_VSITE2, 0, plist);
694 print_bondeds(out,at->nr,d_vsites3, F_VSITE3, 0, plist);
695 print_bondeds(out,at->nr,d_vsites3, F_VSITE3FD, 0, plist);
696 print_bondeds(out,at->nr,d_vsites3, F_VSITE3FAD,0, plist);
697 print_bondeds(out,at->nr,d_vsites3, F_VSITE3OUT,0, plist);
698 print_bondeds(out,at->nr,d_vsites4, F_VSITE4FD, 0, plist);
699 print_bondeds(out,at->nr,d_vsites4, F_VSITE4FDN, 0, plist);
702 print_top_posre(out,pr);
706 static atom_id search_res_atom(const char *type,int resind,
707 int natom,t_atom at[],
708 char ** const *aname,
709 const char *bondtype,gmx_bool bAllowMissing)
713 for(i=0; (i<natom); i++)
714 if (at[i].resind == resind)
715 return search_atom(type,i,natom,at,aname,bondtype,bAllowMissing);
720 static void do_ssbonds(t_params *ps,int natoms,t_atom atom[],char **aname[],
721 int nssbonds,t_ssbond *ssbonds,gmx_bool bAllowMissing)
726 for(i=0; (i<nssbonds); i++) {
727 ri = ssbonds[i].res1;
728 rj = ssbonds[i].res2;
729 ai = search_res_atom(ssbonds[i].a1,ri,natoms,atom,aname,
730 "special bond",bAllowMissing);
731 aj = search_res_atom(ssbonds[i].a2,rj,natoms,atom,aname,
732 "special bond",bAllowMissing);
733 if ((ai == NO_ATID) || (aj == NO_ATID))
734 gmx_fatal(FARGS,"Trying to make impossible special bond (%s-%s)!",
735 ssbonds[i].a1,ssbonds[i].a2);
736 add_param(ps,ai,aj,NULL,NULL);
740 static void at2bonds(t_params *psb, t_hackblock *hb,
741 int natoms, t_atom atom[], char **aname[],
743 real long_bond_dist, real short_bond_dist,
744 gmx_bool bAllowMissing)
748 real dist2, long_bond_dist2, short_bond_dist2;
751 long_bond_dist2 = sqr(long_bond_dist);
752 short_bond_dist2 = sqr(short_bond_dist);
759 fprintf(stderr,"Making bonds...\n");
761 for(resind=0; (resind < nres) && (i<natoms); resind++) {
762 /* add bonds from list of bonded interactions */
763 for(j=0; j < hb[resind].rb[ebtsBONDS].nb; j++) {
764 /* Unfortunately we can not issue errors or warnings
765 * for missing atoms in bonds, as the hydrogens and terminal atoms
766 * have not been added yet.
768 ai=search_atom(hb[resind].rb[ebtsBONDS].b[j].AI,i,natoms,atom,aname,
770 aj=search_atom(hb[resind].rb[ebtsBONDS].b[j].AJ,i,natoms,atom,aname,
772 if (ai != NO_ATID && aj != NO_ATID) {
773 dist2 = distance2(x[ai],x[aj]);
774 if (dist2 > long_bond_dist2 )
776 fprintf(stderr,"Warning: Long Bond (%d-%d = %g nm)\n",
777 ai+1,aj+1,sqrt(dist2));
779 else if (dist2 < short_bond_dist2 )
781 fprintf(stderr,"Warning: Short Bond (%d-%d = %g nm)\n",
782 ai+1,aj+1,sqrt(dist2));
784 add_param(psb,ai,aj,NULL,hb[resind].rb[ebtsBONDS].b[j].s);
787 /* add bonds from list of hacks (each added atom gets a bond) */
788 while( (i<natoms) && (atom[i].resind == resind) ) {
789 for(j=0; j < hb[resind].nhack; j++)
790 if ( ( hb[resind].hack[j].tp > 0 ||
791 hb[resind].hack[j].oname==NULL ) &&
792 strcmp(hb[resind].hack[j].AI,*(aname[i])) == 0 ) {
793 switch(hb[resind].hack[j].tp) {
794 case 9: /* COOH terminus */
795 add_param(psb,i,i+1,NULL,NULL); /* C-O */
796 add_param(psb,i,i+2,NULL,NULL); /* C-OA */
797 add_param(psb,i+2,i+3,NULL,NULL); /* OA-H */
800 for(k=0; (k<hb[resind].hack[j].nr); k++)
801 add_param(psb,i,i+k+1,NULL,NULL);
806 /* we're now at the start of the next residue */
810 static int pcompar(const void *a, const void *b)
821 return strlen(pb->s) - strlen(pa->s);
826 static void clean_bonds(t_params *ps)
832 /* swap atomnumbers in bond if first larger than second: */
833 for(i=0; (i<ps->nr); i++)
834 if ( ps->param[i].AJ < ps->param[i].AI ) {
836 ps->param[i].AI = ps->param[i].AJ;
841 qsort(ps->param,ps->nr,(size_t)sizeof(ps->param[0]),pcompar);
843 /* remove doubles, keep the first one always. */
845 for(i=1; (i<ps->nr); i++) {
846 if ((ps->param[i].AI != ps->param[j-1].AI) ||
847 (ps->param[i].AJ != ps->param[j-1].AJ) ) {
849 cp_param(&(ps->param[j]),&(ps->param[i]));
854 fprintf(stderr,"Number of bonds was %d, now %d\n",ps->nr,j);
858 fprintf(stderr,"No bonds\n");
861 void print_sums(t_atoms *atoms, gmx_bool bSystem)
874 for(i=0; (i<atoms->nr); i++) {
876 qtot+=atoms->atom[i].q;
879 fprintf(stderr,"Total mass%s %.3f a.m.u.\n",where,m);
880 fprintf(stderr,"Total charge%s %.3f e\n",where,qtot);
883 static void check_restp_type(const char *name,int t1,int t2)
887 gmx_fatal(FARGS,"Residues in one molecule have a different '%s' type: %d and %d",name,t1,t2);
891 static void check_restp_types(t_restp *r0,t_restp *r1)
895 check_restp_type("all dihedrals",r0->bAlldih,r1->bAlldih);
896 check_restp_type("nrexcl",r0->nrexcl,r1->nrexcl);
897 check_restp_type("HH14",r0->HH14,r1->HH14);
898 check_restp_type("remove dihedrals",r0->bRemoveDih,r1->bRemoveDih);
900 for(i=0; i<ebtsNR; i++)
902 check_restp_type(btsNames[i],r0->rb[i].type,r1->rb[i].type);
906 void add_atom_to_restp(t_restp *restp,int resnr,int at_start,const t_hack *hack)
910 const char *Hnum="123456";
914 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
916 *restp->atomname[at_start], resnr, restp->resname);
918 strcpy(buf, hack->nname);
919 buf[strlen(buf)+1]='\0';
922 buf[strlen(buf)]='-';
925 restp->natom += hack->nr;
926 srenew(restp->atom, restp->natom);
927 srenew(restp->atomname, restp->natom);
928 srenew(restp->cgnr, restp->natom);
930 for(k=restp->natom-1; k > at_start+hack->nr; k--)
933 restp->atom [k - hack->nr];
935 restp->atomname[k - hack->nr];
937 restp->cgnr [k - hack->nr];
940 for(k=0; k < hack->nr; k++)
942 /* set counter in atomname */
945 buf[strlen(buf)-1] = Hnum[k];
947 snew( restp->atomname[at_start+1+k], 1);
948 restp->atom [at_start+1+k] = *hack->atom;
949 *restp->atomname[at_start+1+k] = strdup(buf);
950 if ( hack->cgnr != NOTSET )
952 restp->cgnr [at_start+1+k] = hack->cgnr;
956 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
961 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
962 int nrtp, t_restp rtp[],
963 int nres, t_resinfo *resinfo,
965 t_hackblock **ntdb, t_hackblock **ctdb,
972 const char *Hnum="123456";
978 /* first the termini */
979 for(i=0; i<nterpairs; i++) {
980 if (rn[i] >= 0 && ntdb[i] != NULL) {
981 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
983 if (rc[i] >= 0 && ctdb[i] != NULL) {
984 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
988 /* then the whole rtp */
989 for(i=0; i < nres; i++) {
990 /* Here we allow a mismatch of one character when looking for the rtp entry.
991 * For such a mismatch there should be only one mismatching name.
992 * This is mainly useful for small molecules such as ions.
993 * Note that this will usually not work for protein, DNA and RNA,
994 * since there the residue names should be listed in residuetypes.dat
995 * and an error will have been generated earlier in the process.
997 key = *resinfo[i].rtp;
998 snew(resinfo[i].rtp,1);
999 *resinfo[i].rtp = search_rtp(key,nrtp,rtp);
1000 res = get_restp(*resinfo[i].rtp,nrtp,rtp);
1001 copy_t_restp(res, &(*restp)[i]);
1003 /* Check that we do not have different bonded types in one molecule */
1004 check_restp_types(&(*restp)[0],&(*restp)[i]);
1007 for(j=0; j<nterpairs && tern==-1; j++) {
1013 for(j=0; j<nterpairs && terc == -1; j++) {
1018 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb,tern>=0,terc>=0);
1020 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1021 (terc >= 0 && ctdb[terc] == NULL))) {
1022 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.");
1024 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1025 (terc >= 0 && ctdb[terc]->nhack == 0))) {
1026 gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends. Select a proper terminal entry.");
1030 /* now perform t_hack's on t_restp's,
1031 i.e. add's and deletes from termini database will be
1032 added to/removed from residue topology
1033 we'll do this on one big dirty loop, so it won't make easy reading! */
1034 for(i=0; i < nres; i++)
1036 for(j=0; j < (*hb)[i].nhack; j++)
1038 if ( (*hb)[i].hack[j].nr )
1040 /* find atom in restp */
1041 for(l=0; l < (*restp)[i].natom; l++)
1042 if ( ( (*hb)[i].hack[j].oname==NULL &&
1043 strcmp((*hb)[i].hack[j].AI, *(*restp)[i].atomname[l])==0 ) ||
1044 ( (*hb)[i].hack[j].oname!=NULL &&
1045 strcmp((*hb)[i].hack[j].oname,*(*restp)[i].atomname[l])==0 ) )
1047 if (l == (*restp)[i].natom)
1049 /* If we are doing an atom rename only, we don't need
1050 * to generate a fatal error if the old name is not found
1053 /* Deleting can happen also only on the input atoms,
1054 * not necessarily always on the rtp entry.
1056 if (!((*hb)[i].hack[j].oname != NULL &&
1057 (*hb)[i].hack[j].nname != NULL) &&
1058 !((*hb)[i].hack[j].oname != NULL &&
1059 (*hb)[i].hack[j].nname == NULL))
1062 "atom %s not found in buiding block %d%s "
1063 "while combining tdb and rtp",
1064 (*hb)[i].hack[j].oname!=NULL ?
1065 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].AI,
1066 i+1,*resinfo[i].rtp);
1071 if ( (*hb)[i].hack[j].oname==NULL ) {
1073 add_atom_to_restp(&(*restp)[i],resinfo[i].nr,l,
1079 if ( (*hb)[i].hack[j].nname==NULL ) {
1080 /* we're deleting */
1082 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1083 *(*restp)[i].atomname[l],
1084 i+1,(*restp)[i].resname);
1085 /* shift the rest */
1086 (*restp)[i].natom--;
1087 for(k=l; k < (*restp)[i].natom; k++) {
1088 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1089 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1090 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1092 /* give back space */
1093 srenew((*restp)[i].atom, (*restp)[i].natom);
1094 srenew((*restp)[i].atomname, (*restp)[i].natom);
1095 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1096 } else { /* nname != NULL */
1097 /* we're replacing */
1099 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1100 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1101 i+1,(*restp)[i].resname);
1102 snew( (*restp)[i].atomname[l], 1);
1103 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1104 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1105 if ( (*hb)[i].hack[j].cgnr != NOTSET )
1106 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1115 static gmx_bool atomname_cmp_nr(const char *anm,t_hack *hack,int *nr)
1122 return (gmx_strcasecmp(anm,hack->nname) == 0);
1126 if (isdigit(anm[strlen(anm)-1]))
1128 *nr = anm[strlen(anm)-1] - '0';
1134 if (*nr <= 0 || *nr > hack->nr)
1140 return (strlen(anm) == strlen(hack->nname) + 1 &&
1141 gmx_strncasecmp(anm,hack->nname,strlen(hack->nname)) == 0);
1146 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba,rvec *x,int atind,
1147 t_restp *rptr,t_hackblock *hbr,
1154 char *start_at,buf[STRLEN];
1156 gmx_bool bReplaceReplace,bFoundInAdd;
1159 oldnm = *pdba->atomname[atind];
1160 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1163 for(j=0; j<hbr->nhack; j++)
1165 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1166 gmx_strcasecmp(oldnm,hbr->hack[j].oname) == 0)
1168 /* This is a replace entry. */
1169 /* Check if we are not replacing a replaced atom. */
1170 bReplaceReplace = FALSE;
1171 for(k=0; k<hbr->nhack; k++) {
1173 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1174 gmx_strcasecmp(hbr->hack[k].nname,hbr->hack[j].oname) == 0)
1176 /* The replace in hack[j] replaces an atom that
1177 * was already replaced in hack[k], we do not want
1178 * second or higher level replaces at this stage.
1180 bReplaceReplace = TRUE;
1183 if (bReplaceReplace)
1185 /* Skip this replace. */
1189 /* This atom still has the old name, rename it */
1190 newnm = hbr->hack[j].nname;
1191 for(k=0; k<rptr->natom; k++)
1193 if (gmx_strcasecmp(newnm,*rptr->atomname[k]) == 0)
1198 if (k == rptr->natom)
1200 /* The new name is not present in the rtp.
1201 * We need to apply the replace also to the rtp entry.
1204 /* We need to find the add hack that can add this atom
1205 * to find out after which atom it should be added.
1207 bFoundInAdd = FALSE;
1208 for(k=0; k<hbr->nhack; k++)
1210 if (hbr->hack[k].oname == NULL &&
1211 hbr->hack[k].nname != NULL &&
1212 atomname_cmp_nr(newnm,&hbr->hack[k],&anmnr))
1216 start_at = hbr->hack[k].a[0];
1220 sprintf(buf,"%s%d",hbr->hack[k].nname,anmnr-1);
1223 for(start_nr=0; start_nr<rptr->natom; start_nr++)
1225 if (gmx_strcasecmp(start_at,(*rptr->atomname[start_nr])) == 0)
1230 if (start_nr == rptr->natom)
1232 gmx_fatal(FARGS,"Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1233 start_at,rptr->resname,newnm);
1235 /* We can add the atom after atom start_nr */
1236 add_atom_to_restp(rptr,resnr,start_nr,
1245 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'",
1247 hbr->hack[j].oname,hbr->hack[j].nname,
1254 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1255 oldnm,rptr->resname,resnr,newnm);
1257 /* Rename the atom in pdba */
1258 snew(pdba->atomname[atind],1);
1259 *pdba->atomname[atind] = strdup(newnm);
1261 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1262 gmx_strcasecmp(oldnm,hbr->hack[j].oname) == 0)
1264 /* This is a delete entry, check if this atom is present
1265 * in the rtp entry of this residue.
1267 for(k=0; k<rptr->natom; k++)
1269 if (gmx_strcasecmp(oldnm,*rptr->atomname[k]) == 0)
1274 if (k == rptr->natom)
1276 /* This atom is not present in the rtp entry,
1277 * delete is now from the input pdba.
1281 printf("Deleting atom '%s' in residue '%s' %d\n",
1282 oldnm,rptr->resname,resnr);
1284 /* We should free the atom name,
1285 * but it might be used multiple times in the symtab.
1286 * sfree(pdba->atomname[atind]);
1288 for(k=atind+1; k<pdba->nr; k++)
1290 pdba->atom[k-1] = pdba->atom[k];
1291 pdba->atomname[k-1] = pdba->atomname[k];
1292 copy_rvec(x[k],x[k-1]);
1303 void match_atomnames_with_rtp(t_restp restp[],t_hackblock hb[],
1304 t_atoms *pdba,rvec *x,
1313 char *start_at,buf[STRLEN];
1315 gmx_bool bFoundInAdd;
1317 for(i=0; i<pdba->nr; i++)
1319 oldnm = *pdba->atomname[i];
1320 resnr = pdba->resinfo[pdba->atom[i].resind].nr;
1321 rptr = &restp[pdba->atom[i].resind];
1322 for(j=0; (j<rptr->natom); j++)
1324 if (gmx_strcasecmp(oldnm,*(rptr->atomname[j])) == 0)
1329 if (j == rptr->natom)
1331 /* Not found yet, check if we have to rename this atom */
1332 if (match_atomnames_with_rtp_atom(pdba,x,i,
1333 rptr,&(hb[pdba->atom[i].resind]),
1336 /* We deleted this atom, decrease the atom counter by 1. */
1343 #define NUM_CMAP_ATOMS 5
1344 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
1348 int natoms = atoms->nr;
1349 t_atom *atom = atoms->atom;
1350 char ***aname = atoms->atomname;
1351 t_resinfo *resinfo = atoms->resinfo;
1352 int nres = atoms->nres;
1354 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1355 int cmap_chainnum=-1, this_residue_index;
1362 fprintf(stderr,"Making cmap torsions...");
1364 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1365 * therefore we get a valgrind invalid 4 byte read error with atom am */
1366 for(residx=0; residx<nres-1; residx++)
1368 /* Add CMAP terms from the list of CMAP interactions */
1369 for(j=0;j<restp[residx].rb[ebtsCMAP].nb; j++)
1372 /* Loop over atoms in a candidate CMAP interaction and
1373 * check that they exist, are from the same chain and are
1374 * from residues labelled as protein. */
1375 for(k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1377 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1378 i,natoms,atom,aname,ptr,TRUE);
1379 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1382 /* This break is necessary, because cmap_atomid[k]
1383 * == NO_ATID cannot be safely used as an index
1384 * into the atom array. */
1387 this_residue_index = atom[cmap_atomid[k]].resind;
1390 cmap_chainnum = resinfo[this_residue_index].chainnum;
1394 /* Does the residue for this atom have the same
1395 * chain number as the residues for previous
1397 bAddCMAP = bAddCMAP &&
1398 cmap_chainnum == resinfo[this_residue_index].chainnum;
1400 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt,*(resinfo[this_residue_index].name));
1405 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);
1411 while(atom[i].resind<residx+1)
1418 /* Start the next residue */
1422 scrub_charge_groups(int *cgnr, int natoms)
1426 for(i=0;i<natoms;i++)
1433 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1434 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1435 int nrtp, t_restp rtp[],
1436 t_restp *restp, t_hackblock *hb,
1437 int nterpairs,t_hackblock **ntdb, t_hackblock **ctdb,
1438 gmx_bool bAllowMissing,
1439 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1440 const char *ff, const char *ffdir,
1442 int nssbonds, t_ssbond *ssbonds,
1443 real long_bond_dist, real short_bond_dist,
1444 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1445 gmx_bool bRenumRes,gmx_bool bRTPresname)
1451 t_params plist[F_NRE];
1458 gmx_residuetype_t rt;
1461 gmx_residuetype_init(&rt);
1464 print_resall(debug, atoms->nres, restp, atype);
1465 dump_hb(debug, atoms->nres, hb);
1469 at2bonds(&(plist[F_BONDS]), hb,
1470 atoms->nr, atoms->atom, atoms->atomname, atoms->nres, *x,
1471 long_bond_dist, short_bond_dist, bAllowMissing);
1473 /* specbonds: disulphide bonds & heme-his */
1474 do_ssbonds(&(plist[F_BONDS]),
1475 atoms->nr, atoms->atom, atoms->atomname, nssbonds, ssbonds,
1478 nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1481 fprintf(stderr,"There were %d missing atoms in molecule %s\n",
1484 gmx_fatal(FARGS,"There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1488 /* Cleanup bonds (sort and rm doubles) */
1489 clean_bonds(&(plist[F_BONDS]));
1491 snew(vsite_type,atoms->nr);
1492 for(i=0; i<atoms->nr; i++)
1493 vsite_type[i]=NOTSET;
1495 /* determine which atoms will be vsites and add dummy masses
1496 also renumber atom numbers in plist[0..F_NRE]! */
1497 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1498 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1501 /* Make Angles and Dihedrals */
1502 fprintf(stderr,"Generating angles, dihedrals and pairs...\n");
1503 snew(excls,atoms->nr);
1504 init_nnb(&nnb,atoms->nr,4);
1505 gen_nnb(&nnb,plist);
1506 print_nnb(&nnb,"NNB");
1507 gen_pad(&nnb,atoms,restp[0].nrexcl,restp[0].HH14,
1508 plist,excls,hb,restp[0].bAlldih,restp[0].bRemoveDih,
1515 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1516 if (plist[F_CMAP].nr > 0)
1518 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1523 /* set mass of all remaining hydrogen atoms */
1525 do_h_mass(&(plist[F_BONDS]),vsite_type,atoms,mHmult,bDeuterate);
1528 /* Cleanup bonds (sort and rm doubles) */
1529 /* clean_bonds(&(plist[F_BONDS]));*/
1532 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1533 " %4d pairs, %4d bonds and %4d virtual sites\n",
1534 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1535 plist[F_LJ14].nr, plist[F_BONDS].nr,
1536 plist[F_VSITE2].nr +
1537 plist[F_VSITE3].nr +
1538 plist[F_VSITE3FD].nr +
1539 plist[F_VSITE3FAD].nr +
1540 plist[F_VSITE3OUT].nr +
1541 plist[F_VSITE4FD].nr +
1542 plist[F_VSITE4FDN].nr );
1544 print_sums(atoms, FALSE);
1546 if (FALSE == bChargeGroups)
1548 scrub_charge_groups(cgnr, atoms->nr);
1553 for(i=0; i<atoms->nres; i++)
1555 atoms->resinfo[i].nr = i + 1;
1556 atoms->resinfo[i].ic = ' ';
1561 fprintf(stderr,"Writing topology\n");
1562 /* We can copy the bonded types from the first restp,
1563 * since the types have to be identical for all residues in one molecule.
1565 for(i=0; i<ebtsNR; i++) {
1566 bts[i] = restp[0].rb[i].type;
1568 write_top(top_file, posre_fn, molname,
1570 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1574 free_t_hackblock(atoms->nres, &hb);
1575 free_t_restp(atoms->nres, &restp);
1576 gmx_residuetype_destroy(rt);
1578 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1580 for (i=0; i<F_NRE; i++)
1581 sfree(plist[i].param);
1582 for (i=0; i<atoms->nr; i++)