2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
63 #include "gpp_nextnb.h"
67 #include "hackblock.h"
69 char atp[7] = "HCNOSX";
70 #define NATP (asize(atp)-1)
72 real blen[NATP][NATP] = {
73 { 0.00, 0.108, 0.105, 0.10, 0.10, 0.10 },
74 { 0.108, 0.15, 0.14, 0.14, 0.16, 0.14 },
75 { 0.105, 0.14, 0.14, 0.14, 0.16, 0.14 },
76 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.14 },
77 { 0.10, 0.16, 0.16, 0.17, 0.20, 0.17 },
78 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.17 }
81 #define MARGIN_FAC 1.1
83 static gmx_bool is_bond(int nnm,t_nm2type nmt[],char *ai,char *aj,real blen)
87 for(i=0; (i<nnm); i++) {
88 for(j=0; (j<nmt[i].nbonds); j++) {
89 if ((((gmx_strncasecmp(ai,nmt[i].elem,1) == 0) &&
90 (gmx_strncasecmp(aj,nmt[i].bond[j],1) == 0)) ||
91 ((gmx_strncasecmp(ai,nmt[i].bond[j],1) == 0) &&
92 (gmx_strncasecmp(aj,nmt[i].elem,1) == 0))) &&
93 (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
100 static int get_atype(char *nm)
104 for(i=0; (i<NATP-1); i++) {
105 if (nm[0] == atp[i]) {
113 void mk_bonds(int nnm,t_nm2type nmt[],
114 t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
115 gmx_bool bPBC,matrix box,gmx_atomprop_t aps)
123 for(i=0; (i<MAXATOMLIST); i++)
125 for(i=0; (i<MAXFORCEPARAM); i++)
129 set_pbc(&pbc,-1,box);
130 for(i=0; (i<atoms->nr); i++) {
132 fprintf(stderr,"\ratom %d",i);
133 for(j=i+1; (j<atoms->nr); j++) {
135 pbc_dx(&pbc,x[i],x[j],dx);
137 rvec_sub(x[i],x[j],dx);
140 if (is_bond(nnm,nmt,*atoms->atomname[i],*atoms->atomname[j],
145 add_param_to_list (bond,&b);
149 fprintf(debug,"Bonding atoms %s-%d and %s-%d\n",
150 *atoms->atomname[i],i+1,*atoms->atomname[j],j+1);
154 fprintf(stderr,"\ratom %d\n",i);
157 int *set_cgnr(t_atoms *atoms,gmx_bool bUsePDBcharge,real *qtot,real *mtot)
164 snew(cgnr,atoms->nr);
165 for(i=0; (i<atoms->nr); i++) {
166 if (atoms->pdbinfo && bUsePDBcharge)
167 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
168 qt += atoms->atom[i].q;
169 *qtot += atoms->atom[i].q;
170 *mtot += atoms->atom[i].m;
180 gpp_atomtype_t set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
181 int *nbonds,int nnm,t_nm2type nm2t[])
183 gpp_atomtype_t atype;
186 atype = init_atomtype();
187 snew(atoms->atomtype,atoms->nr);
188 nresolved = nm2type(nnm,nm2t,tab,atoms,atype,nbonds,bonds);
189 if (nresolved != atoms->nr)
190 gmx_fatal(FARGS,"Could only find a forcefield type for %d out of %d atoms",
191 nresolved,atoms->nr);
193 fprintf(stderr,"There are %d different atom types in your sample\n",
194 get_atomtype_ntypes(atype));
199 void lo_set_force_const(t_params *plist,real c[],int nrfp,gmx_bool bRound,
200 gmx_bool bDih,gmx_bool bParam)
206 for(i=0; (i<plist->nr); i++) {
208 for(j=0; j<nrfp; j++)
212 sprintf(buf,"%.2e",plist->param[i].c[0]);
213 sscanf(buf,"%lf",&cc);
217 c[0] = plist->param[i].c[0];
220 c[0] = ((int)(c[0] + 3600)) % 360;
223 /* To put the minimum at the current angle rather than the maximum */
227 for(j=0; (j<nrfp); j++) {
228 plist->param[i].c[j] = c[j];
229 plist->param[i].c[nrfp+j] = c[j];
231 set_p_string(&(plist->param[i]),"");
235 void set_force_const(t_params plist[],real kb,real kt,real kp,gmx_bool bRound,
239 real c[MAXFORCEPARAM];
243 lo_set_force_const(&plist[F_BONDS],c,2,bRound,FALSE,bParam);
245 lo_set_force_const(&plist[F_ANGLES],c,2,bRound,FALSE,bParam);
248 lo_set_force_const(&plist[F_PDIHS],c,3,bRound,TRUE,bParam);
251 void calc_angles_dihs(t_params *ang,t_params *dih,rvec x[],gmx_bool bPBC,
254 int i,ai,aj,ak,al,t1,t2,t3;
255 rvec r_ij,r_kj,r_kl,m,n;
256 real sign,th,costh,ph;
260 set_pbc(&pbc,epbcXYZ,box);
262 pr_rvecs(debug,0,"X2TOP",box,DIM);
263 for(i=0; (i<ang->nr); i++) {
264 ai = ang->param[i].AI;
265 aj = ang->param[i].AJ;
266 ak = ang->param[i].AK;
267 th = RAD2DEG*bond_angle(x[ai],x[aj],x[ak],bPBC ? &pbc : NULL,
268 r_ij,r_kj,&costh,&t1,&t2);
270 fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
271 ai,aj,ak,norm(r_ij),norm(r_kj),th);
272 ang->param[i].C0 = th;
274 for(i=0; (i<dih->nr); i++) {
275 ai = dih->param[i].AI;
276 aj = dih->param[i].AJ;
277 ak = dih->param[i].AK;
278 al = dih->param[i].AL;
279 ph = RAD2DEG*dih_angle(x[ai],x[aj],x[ak],x[al],bPBC ? & pbc : NULL,
280 r_ij,r_kj,r_kl,m,n,&sign,&t1,&t2,&t3);
282 fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d al=%3d r_ij=%8.3f r_kj=%8.3f r_kl=%8.3f ph=%8.3f\n",
283 ai,aj,ak,al,norm(r_ij),norm(r_kj),norm(r_kl),ph);
284 dih->param[i].C0 = ph;
288 static void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[])
292 for(i=0; (i<atoms->nr); i++) {
293 fprintf(fp,"Atom %5s has %1d bonds\n",*atoms->atomname[i],nbonds[i]);
297 static void print_pl(FILE *fp,t_params plist[],int ftp,const char *name,
302 if (plist[ftp].nr > 0) {
304 fprintf(fp,"[ %s ]\n",name);
305 nral = interaction_function[ftp].nratoms;
306 nrfp = interaction_function[ftp].nrfpA;
307 for(i=0; (i<plist[ftp].nr); i++) {
308 for(j=0; (j<nral); j++)
309 fprintf(fp," %5s",*atomname[plist[ftp].param[i].a[j]]);
310 for(j=0; (j<nrfp); j++) {
311 if (plist[ftp].param[i].c[j] != NOTSET)
312 fprintf(fp," %10.3e",plist[ftp].param[i].c[j]);
319 static void print_rtp(const char *filenm,const char *title,t_atoms *atoms,
320 t_params plist[],gpp_atomtype_t atype,int cgnr[],
327 fp = gmx_fio_fopen(filenm,"w");
328 fprintf(fp,"; %s\n",title);
330 fprintf(fp,"[ %s ]\n",*atoms->resinfo[0].name);
332 fprintf(fp,"[ atoms ]\n");
333 for(i=0; (i<atoms->nr); i++) {
334 tp = atoms->atom[i].type;
335 if ((tpnm = get_atomtype_name(tp,atype)) == NULL)
336 gmx_fatal(FARGS,"tp = %d, i = %d in print_rtp",tp,i);
337 fprintf(fp,"%-8s %12s %8.4f %5d\n",
338 *atoms->atomname[i],tpnm,
339 atoms->atom[i].q,cgnr[i]);
341 print_pl(fp,plist,F_BONDS,"bonds",atoms->atomname);
342 print_pl(fp,plist,F_ANGLES,"angles",atoms->atomname);
343 print_pl(fp,plist,F_PDIHS,"dihedrals",atoms->atomname);
344 print_pl(fp,plist,F_IDIHS,"impropers",atoms->atomname);
349 int cmain(int argc, char *argv[])
351 const char *desc[] = {
352 "[TT]g_x2top[tt] generates a primitive topology from a coordinate file.",
353 "The program assumes all hydrogens are present when defining",
354 "the hybridization from the atom name and the number of bonds.",
355 "The program can also make an [TT].rtp[tt] entry, which you can then add",
356 "to the [TT].rtp[tt] database.[PAR]",
357 "When [TT]-param[tt] is set, equilibrium distances and angles",
358 "and force constants will be printed in the topology for all",
359 "interactions. The equilibrium distances and angles are taken",
360 "from the input coordinates, the force constant are set with",
361 "command line options.",
362 "The force fields somewhat supported currently are:[PAR]",
363 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
364 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
365 "The corresponding data files can be found in the library directory",
366 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
367 "information about file formats. By default, the force field selection",
368 "is interactive, but you can use the [TT]-ff[tt] option to specify",
369 "one of the short names above on the command line instead. In that",
370 "case [TT]g_x2top[tt] just looks for the corresponding file.[PAR]",
372 const char *bugs[] = {
373 "The atom type selection is primitive. Virtually no chemical knowledge is used",
374 "Periodic boundary conditions screw up the bonding",
375 "No improper dihedrals are generated",
376 "The atoms to atomtype translation table is incomplete ([TT]atomname2type.n2t[tt] file in the data directory). Please extend it and send the results back to the GROMACS crew."
379 t_params plist[F_NRE];
381 t_atoms *atoms; /* list with all atoms */
382 gpp_atomtype_t atype;
388 char title[STRLEN],forcefield[32],ffdir[STRLEN];
389 rvec *x; /* coordinates? */
391 int bts[] = { 1,1,1,2 };
392 matrix box; /* box length matrix */
393 int natoms; /* number of atoms in one molecule */
394 int nres; /* number of molecules? */
397 gmx_bool bRTP,bTOP,bOPLS;
399 real cutoff,qtot,mtot;
404 { efSTX, "-f", "conf", ffREAD },
405 { efTOP, "-o", "out", ffOPTWR },
406 { efRTP, "-r", "out", ffOPTWR }
408 #define NFILE asize(fnm)
409 static real scale = 1.1, kb = 4e5,kt = 400,kp = 5;
410 static t_restp rtp_header_settings;
411 static gmx_bool bRemoveDihedralIfWithImproper = FALSE;
412 static gmx_bool bGenerateHH14Interactions = TRUE;
413 static gmx_bool bKeepAllGeneratedDihedrals = FALSE;
414 static int nrexcl = 3;
415 static gmx_bool bParam = TRUE, bRound = TRUE;
416 static gmx_bool bPairs = TRUE, bPBC = TRUE;
417 static gmx_bool bUsePDBcharge = FALSE,bVerbose=FALSE;
418 static const char *molnm = "ICE";
419 static const char *ff = "oplsaa";
421 { "-ff", FALSE, etSTR, {&ff},
422 "Force field for your simulation. Type \"select\" for interactive selection." },
423 { "-v", FALSE, etBOOL, {&bVerbose},
424 "Generate verbose output in the top file." },
425 { "-nexcl", FALSE, etINT, {&nrexcl},
426 "Number of exclusions" },
427 { "-H14", FALSE, etBOOL, {&bGenerateHH14Interactions},
428 "Use 3rd neighbour interactions for hydrogen atoms" },
429 { "-alldih", FALSE, etBOOL, {&bKeepAllGeneratedDihedrals},
430 "Generate all proper dihedrals" },
431 { "-remdih", FALSE, etBOOL, {&bRemoveDihedralIfWithImproper},
432 "Remove dihedrals on the same bond as an improper" },
433 { "-pairs", FALSE, etBOOL, {&bPairs},
434 "Output 1-4 interactions (pairs) in topology file" },
435 { "-name", FALSE, etSTR, {&molnm},
436 "Name of your molecule" },
437 { "-pbc", FALSE, etBOOL, {&bPBC},
438 "Use periodic boundary conditions." },
439 { "-pdbq", FALSE, etBOOL, {&bUsePDBcharge},
440 "Use the B-factor supplied in a [TT].pdb[tt] file for the atomic charges" },
441 { "-param", FALSE, etBOOL, {&bParam},
442 "Print parameters in the output" },
443 { "-round", FALSE, etBOOL, {&bRound},
444 "Round off measured values" },
445 { "-kb", FALSE, etREAL, {&kb},
446 "Bonded force constant (kJ/mol/nm^2)" },
447 { "-kt", FALSE, etREAL, {&kt},
448 "Angle force constant (kJ/mol/rad^2)" },
449 { "-kp", FALSE, etREAL, {&kp},
450 "Dihedral angle force constant (kJ/mol/rad^2)" }
453 CopyRight(stderr,argv[0]);
455 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
456 asize(desc),desc,asize(bugs),bugs,&oenv);
457 bRTP = opt2bSet("-r",NFILE,fnm);
458 bTOP = opt2bSet("-o",NFILE,fnm);
459 /* C89 requirements mean that these struct members cannot be used in
460 * the declaration of pa. So some temporary variables are needed. */
461 rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
462 rtp_header_settings.bGenerateHH14Interactions = bGenerateHH14Interactions;
463 rtp_header_settings.bKeepAllGeneratedDihedrals = bKeepAllGeneratedDihedrals;
464 rtp_header_settings.nrexcl = nrexcl;
467 gmx_fatal(FARGS,"Specify at least one output file");
469 aps = gmx_atomprop_init();
471 /* Force field selection, interactive or direct */
472 choose_ff(strcmp(ff,"select") == 0 ? NULL : ff,
473 forcefield,sizeof(forcefield),
474 ffdir,sizeof(ffdir));
476 bOPLS = (strcmp(forcefield,"oplsaa") == 0);
479 mymol.name = strdup(molnm);
482 /* Init parameter lists */
485 /* Read coordinates */
486 get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
489 /* make space for all the atoms */
490 init_t_atoms(atoms,natoms,TRUE);
493 read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,NULL,&epbc,box);
495 sprintf(n2t,"%s",ffdir);
496 nm2t = rd_nm2type(n2t,&nnm);
498 gmx_fatal(FARGS,"No or incorrect atomname2type.n2t file found (looking for %s)",
501 printf("There are %d name to type translations in file %s\n",nnm,n2t);
503 dump_nm2type(debug,nnm,nm2t);
504 printf("Generating bonds from distances...\n");
505 snew(nbonds,atoms->nr);
506 mk_bonds(nnm,nm2t,atoms,x,&(plist[F_BONDS]),nbonds,forcefield,
509 open_symtab(&symtab);
510 atype = set_atom_type(&symtab,atoms,&(plist[F_BONDS]),nbonds,nnm,nm2t);
512 /* Make Angles and Dihedrals */
513 snew(excls,atoms->nr);
514 printf("Generating angles and dihedrals from bonds...\n");
515 init_nnb(&nnb,atoms->nr,4);
517 print_nnb(&nnb,"NNB");
518 gen_pad(&nnb,atoms,&rtp_header_settings,plist,excls,NULL,TRUE);
522 plist[F_LJ14].nr = 0;
524 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
525 " %4d pairs, %4d bonds and %4d atoms\n",
527 bOPLS ? "Ryckaert-Bellemans" : "proper",
528 plist[F_IDIHS].nr, plist[F_ANGLES].nr,
529 plist[F_LJ14].nr, plist[F_BONDS].nr,atoms->nr);
531 calc_angles_dihs(&plist[F_ANGLES],&plist[F_PDIHS],x,bPBC,box);
533 set_force_const(plist,kb,kt,kp,bRound,bParam);
535 cgnr = set_cgnr(atoms,bUsePDBcharge,&qtot,&mtot);
536 printf("Total charge is %g, total mass is %g\n",qtot,mtot);
543 fp = ftp2FILE(efTOP,NFILE,fnm,"w");
544 print_top_header(fp,ftp2fn(efTOP,NFILE,fnm),
545 "Generated by x2top",TRUE,ffdir,1.0);
547 write_top(fp,NULL,mymol.name,atoms,FALSE,bts,plist,excls,atype,
548 cgnr,rtp_header_settings.nrexcl);
549 print_top_mols(fp,mymol.name,ffdir,NULL,0,NULL,1,&mymol);
554 print_rtp(ftp2fn(efRTP,NFILE,fnm),"Generated by x2top",
555 atoms,plist,atype,cgnr,asize(bts),bts);
558 dump_hybridization(debug,atoms,nbonds);
560 close_symtab(&symtab);
563 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",Program());
564 printf(" Please verify atomtypes and charges by comparison to other\n");
565 printf(" topologies.\n");