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-2004, 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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
49 #include "gmx_fatal.h"
50 #include "gpp_atomtype.h"
54 void set_p_string(t_param
*p
,const char *s
)
57 if (strlen(s
) < sizeof(p
->s
)-1)
58 strncpy(p
->s
,s
,sizeof(p
->s
));
60 gmx_fatal(FARGS
,"Increase MAXSLEN in include/grompp.h to at least %d,"
61 " or shorten your definition of bonds like %s to at most %d",
62 strlen(s
)+1,s
,MAXSLEN
-1);
68 void pr_alloc (int extra
, t_params
*pr
)
72 /* get new space for arrays */
74 gmx_fatal(FARGS
,"Trying to make array smaller.\n");
77 if ((pr
->nr
== 0) && (pr
->param
!= NULL
)) {
78 fprintf(stderr
,"Warning: dangling pointer at %lx\n",
79 (unsigned long)pr
->param
);
82 if (pr
->nr
+extra
> pr
->maxnr
) {
83 pr
->maxnr
= max(1.2*pr
->maxnr
,pr
->maxnr
+ extra
);
84 srenew(pr
->param
,pr
->maxnr
);
85 for(i
=pr
->nr
; (i
<pr
->maxnr
); i
++) {
86 for(j
=0; (j
<MAXATOMLIST
); j
++)
88 for(j
=0; (j
<MAXFORCEPARAM
); j
++)
90 set_p_string(&(pr
->param
[i
]),"");
95 void init_plist(t_params plist
[])
99 for(i
=0; (i
<F_NRE
); i
++) {
102 plist
[i
].param
= NULL
;
107 plist
[i
].grid_spacing
= 0;
110 plist
[i
].cmap_types
= NULL
;
114 void cp_param(t_param
*dest
,t_param
*src
)
118 for(j
=0; (j
<MAXATOMLIST
); j
++)
119 dest
->a
[j
] = src
->a
[j
];
120 for(j
=0; (j
<MAXFORCEPARAM
); j
++)
121 dest
->c
[j
] = src
->c
[j
];
122 strncpy(dest
->s
,src
->s
,sizeof(dest
->s
));
125 void add_param_to_list(t_params
*list
, t_param
*b
)
129 /* allocate one position extra */
132 /* fill the arrays */
133 for (j
=0; (j
< MAXFORCEPARAM
); j
++)
134 list
->param
[list
->nr
].c
[j
] = b
->c
[j
];
135 for (j
=0; (j
< MAXATOMLIST
); j
++)
136 list
->param
[list
->nr
].a
[j
] = b
->a
[j
];
137 memset(list
->param
[list
->nr
].s
,0,sizeof(list
->param
[list
->nr
].s
));
143 void init_molinfo(t_molinfo
*mol
)
146 mol
->excl_set
= FALSE
;
147 mol
->bProcessed
= FALSE
;
148 init_plist(mol
->plist
);
149 init_block(&mol
->cgs
);
150 init_block(&mol
->mols
);
151 init_blocka(&mol
->excls
);
152 init_atom(&mol
->atoms
);
157 void done_bt (t_params
*pl
)
162 void done_mi(t_molinfo
*mi
)
166 done_atom (&(mi
->atoms
));
167 done_block(&(mi
->cgs
));
168 done_block(&(mi
->mols
));
169 for(i
=0; (i
<F_NRE
); i
++)
170 done_bt(&(mi
->plist
[i
]));
173 /* PRINTING STRUCTURES */
175 static void print_nbt (FILE *out
,char *title
,gpp_atomtype_t at
,
176 int ftype
,t_params
*nbt
)
178 int f
,i
,j
,k
,l
,nrfp
,ntype
;
188 fprintf (out
,"; %s\n",title
);
189 fprintf (out
,"[ %s ]\n",dir2str(d_nonbond_params
));
190 fprintf (out
,"; %6s %8s","ai","aj");
191 fprintf (out
," %8s","funct");
192 for (j
=0; (j
<nrfp
); j
++)
193 fprintf (out
," %11c%1d",'c',j
);
196 /* print non-bondtypes */
197 ntype
= get_atomtype_ntypes(at
);
198 for (i
=k
=0; (i
<ntype
); i
++)
199 for(j
=0; (j
<ntype
); j
++,k
++) {
200 fprintf (out
,"%8s %8s %8d",
201 get_atomtype_name(i
,at
),get_atomtype_name(f
,at
),f
);
202 for(l
=0; (l
<nrfp
); l
++)
203 fprintf (out
," %12.5e",nbt
->param
[k
].c
[l
]);
210 void print_bt(FILE *out
, directive d
, gpp_atomtype_t at
,
211 int ftype
,int fsubtype
,t_params plist
[],
214 /* This dihp is a DIRTY patch because the dih-types do not use
215 * all four atoms to determine the type.
217 const int dihp
[2][2] = { { 1,2 }, { 0,3 } };
220 bool bDih
=FALSE
,bSwapParity
;
247 case F_CROSS_BOND_ANGLES
:
250 case F_CROSS_BOND_BONDS
:
296 fprintf(out
,"[ %s ]\n",dir2str(d
));
299 fprintf (out
,"%3s %4s","ai","aj");
300 for (j
=2; (j
<nral
); j
++)
301 fprintf (out
," %3c%c",'a','i'+j
);
304 for (j
=0; (j
<2); j
++)
305 fprintf (out
,"%3c%c",'a','i'+dihp
[f
][j
]);
307 fprintf (out
," funct");
308 for (j
=0; (j
<nrfp
); j
++)
309 fprintf (out
," %12c%1d",'c',j
);
312 /* print bondtypes */
313 for (i
=0; (i
<bt
->nr
); i
++) {
314 bSwapParity
= (bt
->param
[i
].C0
==NOTSET
) && (bt
->param
[i
].C1
==-1);
316 for (j
=0; (j
<nral
); j
++)
317 fprintf (out
,"%5s ",get_atomtype_name(bt
->param
[i
].a
[j
],at
));
320 fprintf (out
,"%5s ",get_atomtype_name(bt
->param
[i
].a
[dihp
[f
][j
]],at
));
321 fprintf (out
,"%5d ", bSwapParity
? -f
-1 : f
+1);
323 if (bt
->param
[i
].s
[0])
324 fprintf(out
," %s",bt
->param
[i
].s
);
326 for (j
=0; (j
<nrfp
&& (bt
->param
[i
].c
[j
] != NOTSET
)); j
++)
327 fprintf (out
,"%13.6e ",bt
->param
[i
].c
[j
]);
335 void print_blocka(FILE *out
, const char *szName
,
336 const char *szIndex
, const char *szA
,
341 fprintf (out
,"; %s\n",szName
);
342 fprintf (out
,"; %4s %s\n",szIndex
,szA
);
343 for (i
=0; (i
< block
->nr
); i
++) {
344 for (i
=0; (i
< block
->nr
); i
++) {
345 fprintf (out
,"%6d",i
+1);
346 for (j
=block
->index
[i
]; (j
< ((int)block
->index
[i
+1])); j
++)
347 fprintf (out
,"%5d",block
->a
[j
]+1);
354 void print_excl(FILE *out
, int natoms
, t_excls excls
[])
361 for(i
=0; i
<natoms
&& !have_excl
; i
++)
362 have_excl
= (excls
[i
].nr
> 0);
365 fprintf (out
,"[ %s ]\n",dir2str(d_exclusions
));
366 fprintf (out
,"; %4s %s\n","i","excluded from i");
367 for(i
=0; i
<natoms
; i
++)
368 if (excls
[i
].nr
> 0) {
369 fprintf (out
,"%6d ",i
+1);
370 for(j
=0; j
<excls
[i
].nr
; j
++)
371 fprintf (out
," %5d",excls
[i
].e
[j
]+1);
379 static double get_residue_charge(const t_atoms
*atoms
,int at
)
384 ri
= atoms
->atom
[at
].resind
;
386 while (at
< atoms
->nr
&& atoms
->atom
[at
].resind
== ri
) {
387 q
+= atoms
->atom
[at
].q
;
394 void print_atoms(FILE *out
,gpp_atomtype_t atype
,t_atoms
*at
,int *cgnr
,
404 fprintf(out
,"[ %s ]\n",as
);
405 fprintf(out
,"; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
406 "nr","type","resnr","residue","atom","cgnr","charge","mass","typeB","chargeB","massB");
411 fprintf(debug
,"This molecule has %d atoms and %d residues\n",
415 /* if the information is present... */
416 for (i
=0; (i
< at
->nr
); i
++) {
417 ri
= at
->atom
[i
].resind
;
418 if ((i
== 0 || ri
!= at
->atom
[i
-1].resind
) &&
419 at
->resinfo
[ri
].rtp
!= NULL
) {
420 qres
= get_residue_charge(at
,i
);
421 fprintf(out
,"; residue %3d %-3s rtp %-4s q ",
423 *at
->resinfo
[ri
].name
,
424 *at
->resinfo
[ri
].rtp
);
425 if (fabs(qres
) < 0.001) {
426 fprintf(out
," %s","0.0");
428 fprintf(out
,"%+3.1f",qres
);
432 tpA
= at
->atom
[i
].type
;
433 if ((tpnmA
= get_atomtype_name(tpA
,atype
)) == NULL
)
434 gmx_fatal(FARGS
,"tpA = %d, i= %d in print_atoms",tpA
,i
);
436 fprintf(out
,"%6d %10s %6d%c %5s %6s %6d %10g %10g",
441 *(at
->resinfo
[at
->atom
[i
].resind
].rtp
) :
442 *(at
->resinfo
[at
->atom
[i
].resind
].name
),
443 *(at
->atomname
[i
]),cgnr
[i
],
444 at
->atom
[i
].q
,at
->atom
[i
].m
);
445 if (PERTURBED(at
->atom
[i
])) {
446 tpB
= at
->atom
[i
].typeB
;
447 if ((tpnmB
= get_atomtype_name(tpB
,atype
)) == NULL
)
448 gmx_fatal(FARGS
,"tpB = %d, i= %d in print_atoms",tpB
,i
);
449 fprintf(out
," %6s %10g %10g",
450 tpnmB
,at
->atom
[i
].qB
,at
->atom
[i
].mB
);
452 qtot
+= (double)at
->atom
[i
].q
;
453 if ( fabs(qtot
) < 4*GMX_REAL_EPS
)
455 fprintf(out
," ; qtot %.4g\n",qtot
);
462 void print_bondeds(FILE *out
,int natoms
,directive d
,
463 int ftype
,int fsubtype
,t_params plist
[])
466 gpp_atomtype_t atype
;
471 atype
= init_atomtype();
475 for (i
=0; (i
< natoms
); i
++) {
477 sprintf(buf
,"%4d",(i
+1));
478 add_atomtype(atype
,&stab
,a
,buf
,param
,0,0,0,0,0,0,0);
480 print_bt(out
,d
,atype
,ftype
,fsubtype
,plist
,TRUE
);
485 done_atomtype(atype
);