3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
60 #include "gpp_nextnb.h"
64 #include "hackblock.h"
66 char atp[7] = "HCNOSX";
67 #define NATP (asize(atp)-1)
69 real blen[NATP][NATP] = {
70 { 0.00, 0.108, 0.105, 0.10, 0.10, 0.10 },
71 { 0.108, 0.15, 0.14, 0.14, 0.16, 0.14 },
72 { 0.105, 0.14, 0.14, 0.14, 0.16, 0.14 },
73 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.14 },
74 { 0.10, 0.16, 0.16, 0.17, 0.20, 0.17 },
75 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.17 }
78 #define MARGIN_FAC 1.1
80 static gmx_bool is_bond(int nnm,t_nm2type nmt[],char *ai,char *aj,real blen)
84 for(i=0; (i<nnm); i++) {
85 for(j=0; (j<nmt[i].nbonds); j++) {
86 if ((((gmx_strncasecmp(ai,nmt[i].elem,1) == 0) &&
87 (gmx_strncasecmp(aj,nmt[i].bond[j],1) == 0)) ||
88 ((gmx_strncasecmp(ai,nmt[i].bond[j],1) == 0) &&
89 (gmx_strncasecmp(aj,nmt[i].elem,1) == 0))) &&
90 (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
97 static int get_atype(char *nm)
101 for(i=0; (i<NATP-1); i++) {
102 if (nm[0] == atp[i]) {
110 void mk_bonds(int nnm,t_nm2type nmt[],
111 t_atoms *atoms,rvec x[],t_params *bond,int nbond[],char *ff,
112 gmx_bool bPBC,matrix box,gmx_atomprop_t aps)
120 for(i=0; (i<MAXATOMLIST); i++)
122 for(i=0; (i<MAXFORCEPARAM); i++)
126 set_pbc(&pbc,-1,box);
127 for(i=0; (i<atoms->nr); i++) {
129 fprintf(stderr,"\ratom %d",i);
130 for(j=i+1; (j<atoms->nr); j++) {
132 pbc_dx(&pbc,x[i],x[j],dx);
134 rvec_sub(x[i],x[j],dx);
137 if (is_bond(nnm,nmt,*atoms->atomname[i],*atoms->atomname[j],
142 add_param_to_list (bond,&b);
146 fprintf(debug,"Bonding atoms %s-%d and %s-%d\n",
147 *atoms->atomname[i],i+1,*atoms->atomname[j],j+1);
151 fprintf(stderr,"\ratom %d\n",i);
154 int *set_cgnr(t_atoms *atoms,gmx_bool bUsePDBcharge,real *qtot,real *mtot)
161 snew(cgnr,atoms->nr);
162 for(i=0; (i<atoms->nr); i++) {
163 if (atoms->pdbinfo && bUsePDBcharge)
164 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
165 qt += atoms->atom[i].q;
166 *qtot += atoms->atom[i].q;
167 *mtot += atoms->atom[i].m;
177 gpp_atomtype_t set_atom_type(t_symtab *tab,t_atoms *atoms,t_params *bonds,
178 int *nbonds,int nnm,t_nm2type nm2t[])
180 gpp_atomtype_t atype;
183 atype = init_atomtype();
184 snew(atoms->atomtype,atoms->nr);
185 nresolved = nm2type(nnm,nm2t,tab,atoms,atype,nbonds,bonds);
186 if (nresolved != atoms->nr)
187 gmx_fatal(FARGS,"Could only find a forcefield type for %d out of %d atoms",
188 nresolved,atoms->nr);
190 fprintf(stderr,"There are %d different atom types in your sample\n",
191 get_atomtype_ntypes(atype));
196 void lo_set_force_const(t_params *plist,real c[],int nrfp,gmx_bool bRound,
197 gmx_bool bDih,gmx_bool bParam)
203 for(i=0; (i<plist->nr); i++) {
205 for(j=0; j<nrfp; j++)
209 sprintf(buf,"%.2e",plist->param[i].c[0]);
210 sscanf(buf,"%lf",&cc);
214 c[0] = plist->param[i].c[0];
217 c[0] = ((int)(c[0] + 3600)) % 360;
220 /* To put the minimum at the current angle rather than the maximum */
224 for(j=0; (j<nrfp); j++) {
225 plist->param[i].c[j] = c[j];
226 plist->param[i].c[nrfp+j] = c[j];
228 set_p_string(&(plist->param[i]),"");
232 void set_force_const(t_params plist[],real kb,real kt,real kp,gmx_bool bRound,
236 real c[MAXFORCEPARAM];
240 lo_set_force_const(&plist[F_BONDS],c,2,bRound,FALSE,bParam);
242 lo_set_force_const(&plist[F_ANGLES],c,2,bRound,FALSE,bParam);
245 lo_set_force_const(&plist[F_PDIHS],c,3,bRound,TRUE,bParam);
248 void calc_angles_dihs(t_params *ang,t_params *dih,rvec x[],gmx_bool bPBC,
251 int i,ai,aj,ak,al,t1,t2,t3;
252 rvec r_ij,r_kj,r_kl,m,n;
253 real sign,th,costh,ph;
257 set_pbc(&pbc,epbcXYZ,box);
259 pr_rvecs(debug,0,"X2TOP",box,DIM);
260 for(i=0; (i<ang->nr); i++) {
261 ai = ang->param[i].AI;
262 aj = ang->param[i].AJ;
263 ak = ang->param[i].AK;
264 th = RAD2DEG*bond_angle(x[ai],x[aj],x[ak],bPBC ? &pbc : NULL,
265 r_ij,r_kj,&costh,&t1,&t2);
267 fprintf(debug,"X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
268 ai,aj,ak,norm(r_ij),norm(r_kj),th);
269 ang->param[i].C0 = th;
271 for(i=0; (i<dih->nr); i++) {
272 ai = dih->param[i].AI;
273 aj = dih->param[i].AJ;
274 ak = dih->param[i].AK;
275 al = dih->param[i].AL;
276 ph = RAD2DEG*dih_angle(x[ai],x[aj],x[ak],x[al],bPBC ? & pbc : NULL,
277 r_ij,r_kj,r_kl,m,n,&sign,&t1,&t2,&t3);
279 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",
280 ai,aj,ak,al,norm(r_ij),norm(r_kj),norm(r_kl),ph);
281 dih->param[i].C0 = ph;
285 static void dump_hybridization(FILE *fp,t_atoms *atoms,int nbonds[])
289 for(i=0; (i<atoms->nr); i++) {
290 fprintf(fp,"Atom %5s has %1d bonds\n",*atoms->atomname[i],nbonds[i]);
294 static void print_pl(FILE *fp,t_params plist[],int ftp,const char *name,
299 if (plist[ftp].nr > 0) {
301 fprintf(fp,"[ %s ]\n",name);
302 nral = interaction_function[ftp].nratoms;
303 nrfp = interaction_function[ftp].nrfpA;
304 for(i=0; (i<plist[ftp].nr); i++) {
305 for(j=0; (j<nral); j++)
306 fprintf(fp," %5s",*atomname[plist[ftp].param[i].a[j]]);
307 for(j=0; (j<nrfp); j++) {
308 if (plist[ftp].param[i].c[j] != NOTSET)
309 fprintf(fp," %10.3e",plist[ftp].param[i].c[j]);
316 static void print_rtp(const char *filenm,const char *title,t_atoms *atoms,
317 t_params plist[],gpp_atomtype_t atype,int cgnr[],
324 fp = gmx_fio_fopen(filenm,"w");
325 fprintf(fp,"; %s\n",title);
327 fprintf(fp,"[ %s ]\n",*atoms->resinfo[0].name);
329 fprintf(fp,"[ atoms ]\n");
330 for(i=0; (i<atoms->nr); i++) {
331 tp = atoms->atom[i].type;
332 if ((tpnm = get_atomtype_name(tp,atype)) == NULL)
333 gmx_fatal(FARGS,"tp = %d, i = %d in print_rtp",tp,i);
334 fprintf(fp,"%-8s %12s %8.4f %5d\n",
335 *atoms->atomname[i],tpnm,
336 atoms->atom[i].q,cgnr[i]);
338 print_pl(fp,plist,F_BONDS,"bonds",atoms->atomname);
339 print_pl(fp,plist,F_ANGLES,"angles",atoms->atomname);
340 print_pl(fp,plist,F_PDIHS,"dihedrals",atoms->atomname);
341 print_pl(fp,plist,F_IDIHS,"impropers",atoms->atomname);
346 int cmain(int argc, char *argv[])
348 const char *desc[] = {
349 "[TT]g_x2top[tt] generates a primitive topology from a coordinate file.",
350 "The program assumes all hydrogens are present when defining",
351 "the hybridization from the atom name and the number of bonds.",
352 "The program can also make an [TT].rtp[tt] entry, which you can then add",
353 "to the [TT].rtp[tt] database.[PAR]",
354 "When [TT]-param[tt] is set, equilibrium distances and angles",
355 "and force constants will be printed in the topology for all",
356 "interactions. The equilibrium distances and angles are taken",
357 "from the input coordinates, the force constant are set with",
358 "command line options.",
359 "The force fields somewhat supported currently are:[PAR]",
360 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
361 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
362 "The corresponding data files can be found in the library directory",
363 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
364 "information about file formats. By default, the force field selection",
365 "is interactive, but you can use the [TT]-ff[tt] option to specify",
366 "one of the short names above on the command line instead. In that",
367 "case [TT]g_x2top[tt] just looks for the corresponding file.[PAR]",
369 const char *bugs[] = {
370 "The atom type selection is primitive. Virtually no chemical knowledge is used",
371 "Periodic boundary conditions screw up the bonding",
372 "No improper dihedrals are generated",
373 "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."
376 t_params plist[F_NRE];
378 t_atoms *atoms; /* list with all atoms */
379 gpp_atomtype_t atype;
385 char title[STRLEN],forcefield[32],ffdir[STRLEN];
386 rvec *x; /* coordinates? */
388 int bts[] = { 1,1,1,2 };
389 matrix box; /* box length matrix */
390 int natoms; /* number of atoms in one molecule */
391 int nres; /* number of molecules? */
394 gmx_bool bRTP,bTOP,bOPLS;
396 real cutoff,qtot,mtot;
401 { efSTX, "-f", "conf", ffREAD },
402 { efTOP, "-o", "out", ffOPTWR },
403 { efRTP, "-r", "out", ffOPTWR }
405 #define NFILE asize(fnm)
406 static real scale = 1.1, kb = 4e5,kt = 400,kp = 5;
407 static t_restp rtp_header_settings;
408 static gmx_bool bRemoveDihedralIfWithImproper = FALSE;
409 static gmx_bool bGenerateHH14Interactions = TRUE;
410 static gmx_bool bKeepAllGeneratedDihedrals = FALSE;
411 static int nrexcl = 3;
412 static gmx_bool bParam = TRUE, bRound = TRUE;
413 static gmx_bool bPairs = TRUE, bPBC = TRUE;
414 static gmx_bool bUsePDBcharge = FALSE,bVerbose=FALSE;
415 static const char *molnm = "ICE";
416 static const char *ff = "oplsaa";
418 { "-ff", FALSE, etSTR, {&ff},
419 "Force field for your simulation. Type \"select\" for interactive selection." },
420 { "-v", FALSE, etBOOL, {&bVerbose},
421 "Generate verbose output in the top file." },
422 { "-nexcl", FALSE, etINT, {&nrexcl},
423 "Number of exclusions" },
424 { "-H14", FALSE, etBOOL, {&bGenerateHH14Interactions},
425 "Use 3rd neighbour interactions for hydrogen atoms" },
426 { "-alldih", FALSE, etBOOL, {&bKeepAllGeneratedDihedrals},
427 "Generate all proper dihedrals" },
428 { "-remdih", FALSE, etBOOL, {&bRemoveDihedralIfWithImproper},
429 "Remove dihedrals on the same bond as an improper" },
430 { "-pairs", FALSE, etBOOL, {&bPairs},
431 "Output 1-4 interactions (pairs) in topology file" },
432 { "-name", FALSE, etSTR, {&molnm},
433 "Name of your molecule" },
434 { "-pbc", FALSE, etBOOL, {&bPBC},
435 "Use periodic boundary conditions." },
436 { "-pdbq", FALSE, etBOOL, {&bUsePDBcharge},
437 "Use the B-factor supplied in a [TT].pdb[tt] file for the atomic charges" },
438 { "-param", FALSE, etBOOL, {&bParam},
439 "Print parameters in the output" },
440 { "-round", FALSE, etBOOL, {&bRound},
441 "Round off measured values" },
442 { "-kb", FALSE, etREAL, {&kb},
443 "Bonded force constant (kJ/mol/nm^2)" },
444 { "-kt", FALSE, etREAL, {&kt},
445 "Angle force constant (kJ/mol/rad^2)" },
446 { "-kp", FALSE, etREAL, {&kp},
447 "Dihedral angle force constant (kJ/mol/rad^2)" }
450 CopyRight(stderr,argv[0]);
452 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
453 asize(desc),desc,asize(bugs),bugs,&oenv);
454 bRTP = opt2bSet("-r",NFILE,fnm);
455 bTOP = opt2bSet("-o",NFILE,fnm);
456 /* C89 requirements mean that these struct members cannot be used in
457 * the declaration of pa. So some temporary variables are needed. */
458 rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
459 rtp_header_settings.bGenerateHH14Interactions = bGenerateHH14Interactions;
460 rtp_header_settings.bKeepAllGeneratedDihedrals = bKeepAllGeneratedDihedrals;
461 rtp_header_settings.nrexcl = nrexcl;
464 gmx_fatal(FARGS,"Specify at least one output file");
466 aps = gmx_atomprop_init();
468 /* Force field selection, interactive or direct */
469 choose_ff(strcmp(ff,"select") == 0 ? NULL : ff,
470 forcefield,sizeof(forcefield),
471 ffdir,sizeof(ffdir));
473 bOPLS = (strcmp(forcefield,"oplsaa") == 0);
476 mymol.name = strdup(molnm);
479 /* Init parameter lists */
482 /* Read coordinates */
483 get_stx_coordnum(opt2fn("-f",NFILE,fnm),&natoms);
486 /* make space for all the atoms */
487 init_t_atoms(atoms,natoms,TRUE);
490 read_stx_conf(opt2fn("-f",NFILE,fnm),title,atoms,x,NULL,&epbc,box);
492 sprintf(n2t,"%s",ffdir);
493 nm2t = rd_nm2type(n2t,&nnm);
495 gmx_fatal(FARGS,"No or incorrect atomname2type.n2t file found (looking for %s)",
498 printf("There are %d name to type translations in file %s\n",nnm,n2t);
500 dump_nm2type(debug,nnm,nm2t);
501 printf("Generating bonds from distances...\n");
502 snew(nbonds,atoms->nr);
503 mk_bonds(nnm,nm2t,atoms,x,&(plist[F_BONDS]),nbonds,forcefield,
506 open_symtab(&symtab);
507 atype = set_atom_type(&symtab,atoms,&(plist[F_BONDS]),nbonds,nnm,nm2t);
509 /* Make Angles and Dihedrals */
510 snew(excls,atoms->nr);
511 printf("Generating angles and dihedrals from bonds...\n");
512 init_nnb(&nnb,atoms->nr,4);
514 print_nnb(&nnb,"NNB");
515 gen_pad(&nnb,atoms,&rtp_header_settings,plist,excls,NULL,TRUE);
519 plist[F_LJ14].nr = 0;
521 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
522 " %4d pairs, %4d bonds and %4d atoms\n",
524 bOPLS ? "Ryckaert-Bellemans" : "proper",
525 plist[F_IDIHS].nr, plist[F_ANGLES].nr,
526 plist[F_LJ14].nr, plist[F_BONDS].nr,atoms->nr);
528 calc_angles_dihs(&plist[F_ANGLES],&plist[F_PDIHS],x,bPBC,box);
530 set_force_const(plist,kb,kt,kp,bRound,bParam);
532 cgnr = set_cgnr(atoms,bUsePDBcharge,&qtot,&mtot);
533 printf("Total charge is %g, total mass is %g\n",qtot,mtot);
540 fp = ftp2FILE(efTOP,NFILE,fnm,"w");
541 print_top_header(fp,ftp2fn(efTOP,NFILE,fnm),
542 "Generated by x2top",TRUE,ffdir,1.0);
544 write_top(fp,NULL,mymol.name,atoms,FALSE,bts,plist,excls,atype,
545 cgnr,rtp_header_settings.nrexcl);
546 print_top_mols(fp,mymol.name,ffdir,NULL,0,NULL,1,&mymol);
551 print_rtp(ftp2fn(efRTP,NFILE,fnm),"Generated by x2top",
552 atoms,plist,atype,cgnr,asize(bts),bts);
555 dump_hybridization(debug,atoms,nbonds);
557 close_symtab(&symtab);
560 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",Program());
561 printf(" Please verify atomtypes and charges by comparison to other\n");
562 printf(" topologies.\n");