← Back to branch summary

~ubuntu-branches/ubuntu/vivid/gdis/vivid

« back to all changes in this revision

Viewing changes to file_gromacs.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Leidert (dale)
  • Date: 2009-04-06 17:12:18 UTC
  • mfrom: (3.1.3 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090406171218-uizoe126jrq09ytt
Tags: 0.90-1
http://bugs.debian.org/492994
https://launchpad.net/bugs/111353
* New upstream release 0.90.
* Acknowledge NMU (closes: #492994). Thanks to Neil Williams.

* makefile.debian: Upstream doesn't provide a Makefile to edit - so created
  our own.
* debian/compat: Added and set to be 5.
* debian/control: Added Homepage, Vcs* and DM-Upload-Allowed fields.
  (Maintainer): Set to the debichem team with approval by Noèl.
  (Uploaders): Added myself.
  (Build-Depends): Increased debhelper version to 5. Removed glutg3-dev.
  Added dpatch.
  (Standards-Version): Bumped to 3.8.1.
  (Description): Removed homepage. Fixed a typo.
* debian/copyright: Updated, completed and adjusted.
* debian/dirs: Dropped useless file.
* debian/docs: Renamed to debian/gdis.docs.
* debian/menu: Renamed to debian/gdis.menu.
  (section): Fixed accordingly to policy.
* debian/gdis.1: Just some formatting changes.
* debian/gdis.desktop: Added file (with small fixes) provided by Phill Bull
  (LP: #111353).
* debian/gdis.install: Added.
* debian/rules: Cleaned. Installation is now done by dh_install. Make sure,
  the CVS directory is not copied. Added dh_desktop call.
* debian/source.lintian-overrides: makefile.debian is created for Debian but
  lives outside debian/.
* debian/watch: Added.
* debian/README.build: Dropped.
* debian/README.source: Added to be compliant to the policy v3.8.
* debian/patches/Debian_make.dpatch: Added.
  - gdis.h (ELEM_FILE): Moved fix for gdis.elemts path (#399132) to this
    patch.
* debian/patches/00list: Added.

Show diffs side-by-side

added added

removed removed

Lines of Context:
file_gromacs.c
 
1
/*
 
2
Copyright (C) 2003 by Sean David Fleming
 
3
 
 
4
sean@ivec.org
 
5
 
 
6
This program is free software; you can redistribute it and/or
 
7
modify it under the terms of the GNU General Public License
 
8
as published by the Free Software Foundation; either version 2
 
9
of the License, or (at your option) any later version.
 
10
 
 
11
This program is distributed in the hope that it will be useful,
 
12
but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
GNU General Public License for more details.
 
15
 
 
16
You should have received a copy of the GNU General Public License
 
17
along with this program; if not, write to the Free Software
 
18
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 
19
 
 
20
The GNU GPL can also be found at http://www.gnu.org
 
21
*/
 
22
 
 
23
#include <math.h>
 
24
#include <stdio.h>
 
25
#include <string.h>
 
26
 
 
27
#include "gdis.h"
 
28
#include "coords.h"
 
29
#include "matrix.h"
 
30
#include "measure.h"
 
31
#include "model.h"
 
32
#include "parse.h"
 
33
#include "file.h"
 
34
#include "render.h"
 
35
#include "scan.h"
 
36
#include "select.h"
 
37
#include "ff.h"
 
38
#include "ff_pak.h"
 
39
#include "interface.h"
 
40
 
 
41
/* main structures */
 
42
extern struct sysenv_pak sysenv;
 
43
extern struct elem_pak elements[];
 
44
 
 
45
/************************************************/
 
46
/* convert internal forcefield type to GROMACS  */
 
47
/************************************************/
 
48
gpointer gromacs_type(struct forcefield_pak *ff)
 
49
{
 
50
struct forcefield_pak *gff;
 
51
 
 
52
gff = ff_dup(ff);
 
53
 
 
54
switch (ff->type)
 
55
  {
 
56
/* bonds */
 
57
  case FF_HARMONIC:
 
58
    switch (ff->data_expected)
 
59
       {
 
60
/* harmonic */
 
61
       case 1:
 
62
         gff->type = 1;
 
63
         gff->data_expected = 1;
 
64
         break;
 
65
/* cubic */
 
66
       case 2:
 
67
       case 3:
 
68
         gff->type = 4;
 
69
         gff->data_expected = 2;
 
70
         break;
 
71
 
 
72
       default:
 
73
         g_assert_not_reached();
 
74
       }
 
75
    break;
 
76
 
 
77
/* morse */
 
78
  case FF_MORSE:
 
79
    gff->type = 2;
 
80
    gff->data_expected = 2;
 
81
    break;
 
82
 
 
83
/* three body */
 
84
  case FF_3B_HARMONIC:
 
85
    switch (ff->data_expected)
 
86
      {
 
87
      case 3:
 
88
        gff->type = 3;
 
89
        gff->data_expected = 5;
 
90
        break;
 
91
 
 
92
      case 1:
 
93
      case 2:
 
94
      default:
 
95
        gff->type = 1;
 
96
        gff->data_expected = 1;
 
97
       }
 
98
    break;
 
99
 
 
100
/* four body */
 
101
  case FF_DIHEDRAL:
 
102
    gff->type = 1;
 
103
    gff->data_expected = 1;
 
104
    break;
 
105
  case FF_IMPROPER:
 
106
    gff->type = 2;
 
107
    gff->data_expected = 1;
 
108
    break;
 
109
  case FF_DIHEDRAL_RB:
 
110
    gff->type = 3;
 
111
    gff->data_expected = 6;
 
112
    break;
 
113
 
 
114
/* non bond */
 
115
  case FF_LENNARD:
 
116
    gff->type = 1;
 
117
    gff->data_expected = 2;
 
118
    break;
 
119
  case FF_BUCKINGHAM:
 
120
    gff->type = 2;
 
121
    gff->data_expected = 2;
 
122
    break;
 
123
  }
 
124
 
 
125
/* GROMACS length units are nm */
 
126
switch (ff->bond_units)
 
127
  {
 
128
  case FF_ANG:
 
129
    gff->bond_value *= 0.1;
 
130
    break;
 
131
  }
 
132
 
 
133
return(gff);
 
134
}
 
135
 
 
136
/******************************/
 
137
/* raw coordinate data (.gro) */
 
138
/******************************/
 
139
gint read_gromacs_gro(gchar *filename, struct model_pak *model)
 
140
{
 
141
gint num_tokens;
 
142
gchar *line, **buff;
 
143
FILE *fp;
 
144
 
 
145
g_assert(model != NULL);
 
146
 
 
147
fp = fopen(filename,"rt");
 
148
if (!fp)
 
149
  return(1);
 
150
 
 
151
/* skip title */
 
152
line = file_read_line(fp);
 
153
g_free(line);
 
154
line = file_read_line(fp);
 
155
 
 
156
while (line)
 
157
  {
 
158
 
 
159
/* FIXME - properly should do formatted read, since things can be */
 
160
/* adjacent with no spaces in between -> tokenize() will fail */
 
161
/*
 
162
  fprintf(fp, "%5i%5s%5s%5i%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n",
 
163
               n, "UNK", core->atom_label, 1, 
 
164
               x[0], x[1], x[2], 0.1*core->v[0], 0.1*core->v[1], 0.1*core->v[2]);
 
165
*/
 
166
 
 
167
  if (strlen(line) > 10)
 
168
    {
 
169
/* tokenize from atom label (avoids some problems with no space separation) */
 
170
    buff = tokenize(&line[10], &num_tokens);
 
171
 
 
172
    if (num_tokens > 4)
 
173
      {
 
174
/* TODO - deal with water properly (ie set atom_types) */
 
175
      if (g_ascii_strncasecmp(*(buff+0), "HW", 2) == 0)
 
176
        {
 
177
        g_free(*(buff+0));
 
178
        *(buff+0) = g_strdup("H");
 
179
        }
 
180
      if (g_ascii_strncasecmp(*(buff+0), "OW", 2) == 0)
 
181
        {
 
182
        g_free(*(buff+0));
 
183
        *(buff+0) = g_strdup("O");
 
184
        }
 
185
 
 
186
      if (elem_test(*(buff+0)))
 
187
        {
 
188
        struct core_pak *core = new_core(*(buff+0), model);
 
189
 
 
190
        core->x[0] = 10.0*str_to_float(*(buff+2));
 
191
        core->x[1] = 10.0*str_to_float(*(buff+3));
 
192
        core->x[2] = 10.0*str_to_float(*(buff+4));
 
193
 
 
194
        model->cores = g_slist_prepend(model->cores, core);
 
195
        }
 
196
      }
 
197
    g_strfreev(buff);
 
198
    }
 
199
  g_free(line);
 
200
  line = file_read_line(fp);
 
201
  }
 
202
 
 
203
model->cores = g_slist_reverse(model->cores);
 
204
 
 
205
model_prep(model);
 
206
 
 
207
return(0);
 
208
}
 
209
 
 
210
/******************************/
 
211
/* raw coordinate data (.gro) */
 
212
/******************************/
 
213
gint write_gromacs_gro(gchar *filename, struct model_pak *model)
 
214
{
 
215
gint n, c, i, s;
 
216
gdouble x[3], min[3], max[3];
 
217
GSList *list;
 
218
struct core_pak *core;
 
219
struct shel_pak *shell;
 
220
FILE *fp;
 
221
 
 
222
fp = fopen(filename, "wt");
 
223
if (!fp)
 
224
  return(1);
 
225
 
 
226
/* init min/max for box calculation */
 
227
for (i=3 ; i-- ; )
 
228
  {
 
229
  max[i] = -G_MAXDOUBLE;
 
230
  min[i] = G_MAXDOUBLE;
 
231
  }
 
232
 
 
233
fprintf(fp, "GROMACS coordinates (generated by GDIS v%f)\n", VERSION);
 
234
 
 
235
c = g_slist_length(model->cores);
 
236
s = g_slist_length(model->shels);
 
237
 
 
238
fprintf(fp, "%5d\n", c+s);
 
239
 
 
240
n = 1;
 
241
for (list=model->cores ; list ; list=g_slist_next(list))
 
242
  {
 
243
  core = list->data;
 
244
 
 
245
/* cartesians */
 
246
  ARR3SET(x, core->x);
 
247
  vecmat(model->latmat, x);
 
248
/* angs -> nm */
 
249
  VEC3MUL(x, 0.1);
 
250
  fprintf(fp, "%5i%5s%5s%5i%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n",
 
251
               n, "UNK", core->atom_label, 1, 
 
252
               x[0], x[1], x[2], 0.1*core->v[0], 0.1*core->v[1], 0.1*core->v[2]);
 
253
 
 
254
/* record min/max for box calculation */
 
255
  for (i=3 ; i-- ; )
 
256
    {
 
257
    if (x[i] < min[i])
 
258
      min[i] = x[i];
 
259
    if (x[i] > max[i])
 
260
      max[i] = x[i];
 
261
    }
 
262
 
 
263
  if (core->shell)
 
264
    {
 
265
/* TODO - has an 'x' post-fix to the atom_label */
 
266
    }
 
267
 
 
268
  n++;
 
269
  }
 
270
 
 
271
/* this is the distance between the min and max atom coords in x, y, and z */
 
272
/* box calculation */
 
273
if (model->periodic == 3)
 
274
  {
 
275
  fprintf(fp, "% f %f %f\n", 0.1*model->pbc[0], 0.1*model->pbc[1], 0.1*model->pbc[2]);
 
276
  }
 
277
else
 
278
  {
 
279
  for (i=0 ; i<3 ; i++)
 
280
    fprintf(fp, " %f", fabs(max[i] - min[i]));
 
281
  fprintf(fp, "\n");
 
282
  }
 
283
 
 
284
fclose(fp);
 
285
 
 
286
return(0);
 
287
}
 
288
 
 
289
/********************/
 
290
/* write the header */
 
291
/********************/
 
292
gint write_gromacs_header(FILE *fp, struct model_pak *model)
 
293
{
 
294
 
 
295
fprintf(fp, "\n;\n; Generated by GDIS v%f\n;\n", VERSION);
 
296
 
 
297
/* put in footer? */
 
298
fprintf(fp, "\n[ system ]\n");
 
299
fprintf(fp, "%s\n", model->basename);
 
300
 
 
301
fprintf(fp, "\n[ moleculetype ]\n");
 
302
/* label, exclusion type */
 
303
fprintf(fp, "%s %d\n", model->basename, 3);
 
304
 
 
305
return(0);
 
306
}
 
307
 
 
308
/********************/
 
309
/* write the footer */
 
310
/********************/
 
311
gint write_gromacs_footer(FILE *fp, struct model_pak *model)
 
312
{
 
313
 
 
314
fprintf(fp, "\n[ molecules ]\n");
 
315
/* label, number of occurences */
 
316
fprintf(fp, "%s 1\n", model->basename);
 
317
 
 
318
return(0);
 
319
}
 
320
 
 
321
/*******************/
 
322
/* write the atoms */
 
323
/*******************/
 
324
gint write_gromacs_atoms(FILE *fp, struct model_pak *model)
 
325
{
 
326
gint n;
 
327
gdouble w;
 
328
GSList *list;
 
329
struct core_pak *core;
 
330
 
 
331
fprintf(fp, "\n[ atoms ]\n");
 
332
 
 
333
n = 1;
 
334
for (list=model->cores ; list ; list=g_slist_next(list))
 
335
  {
 
336
  core = list->data;
 
337
 
 
338
  w = elements[core->atom_code].weight;
 
339
 
 
340
/* if 2nd column is atom type - use if available */
 
341
  if (core->atom_type)
 
342
    {
 
343
    fprintf(fp, "%9d %5s %5d %5s %5s %5d %8.4f %8.4f",
 
344
                 n, core->atom_type, 1, "UNK", core->atom_label, n,
 
345
                 core->charge, w);
 
346
    }
 
347
  else
 
348
    {
 
349
    fprintf(fp, "%9d %5s %5d %5s %5s %5d %8.4f %8.4f",
 
350
                 n, core->atom_label, 1, "UNK", core->atom_label, n,
 
351
                 core->charge, w);
 
352
    }
 
353
 
 
354
/* the lambda = 1.0 values - default to dummy for the time being */
 
355
  fprintf(fp, "    DUM  0.0  %8.4f\n", w);
 
356
 
 
357
  n++;
 
358
  }
 
359
 
 
360
return(0);
 
361
}
 
362
 
 
363
/********************************/
 
364
/* write forcefield information */
 
365
/********************************/
 
366
gint write_gromacs_bonds(FILE *fp, struct model_pak *model)
 
367
{
 
368
gint i, j, n;
 
369
GSList *list;
 
370
struct core_pak *core[2];
 
371
struct forcefield_pak *ff, *ff_gromacs;
 
372
 
 
373
g_assert(model != NULL);
 
374
g_assert(fp != NULL);
 
375
 
 
376
/* header */
 
377
fprintf(fp, "\n[ bonds ]\n");
 
378
 
 
379
/* for each bond FF type in the list -> ff_bond_search -> output */
 
380
for (list=model->bonds ; list ; list=g_slist_next(list))
 
381
  {
 
382
  struct bond_pak *bond = list->data;
 
383
 
 
384
  core[0] = bond->atom1;
 
385
  core[1] = bond->atom2;
 
386
 
 
387
  ff = ff_search(core, 2, model->ff_list);
 
388
 
 
389
  if (ff)
 
390
    {
 
391
/* get core numbering order */
 
392
    i = g_slist_index(model->cores, core[0]) + 1;
 
393
    j = g_slist_index(model->cores, core[1]) + 1;
 
394
 
 
395
    ff_gromacs = gromacs_type(ff);
 
396
 
 
397
/* lambda 0 */
 
398
/* print atoms and type */
 
399
    fprintf(fp, "%2d %2d   %d  ", i, j, ff_gromacs->type);
 
400
 
 
401
/* NB: gromacs has the bond length 1st, followed by the potential parameters */
 
402
/* which is the reverse from internal/GULP format */
 
403
    if (ff_gromacs->bond_expected)
 
404
      fprintf(fp, "%f", ff_gromacs->bond_value);
 
405
    for (n=0 ; n<ff_gromacs->data_expected ; n++)
 
406
      fprintf(fp, " %f", ff_gromacs->data[n]);
 
407
 
 
408
/* lambda 1 */
 
409
    fprintf(fp, "    ");
 
410
    if (ff_gromacs->bond_expected)
 
411
      fprintf(fp, "%f", ff_gromacs->bond_value);
 
412
    for (n=0 ; n<ff_gromacs->data_expected ; n++)
 
413
      fprintf(fp, " %f", ff_gromacs->data[n]);
 
414
 
 
415
    fprintf(fp, "\n");
 
416
 
 
417
    g_free(ff_gromacs);
 
418
    }
 
419
  } 
 
420
 
 
421
return(0);
 
422
}
 
423
 
 
424
/********************************/
 
425
/* write forcefield information */
 
426
/********************************/
 
427
#define MAX_CORES 10
 
428
gint write_gromacs_angles(FILE *fp, struct model_pak *model)
 
429
{
 
430
gint i, j, k, m, n, ni, nj;
 
431
GSList *list1, *list2, *nlist;
 
432
struct core_pak *core[3], *nc[MAX_CORES];
 
433
struct forcefield_pak *ff, *ff_gromacs;
 
434
 
 
435
g_assert(model != NULL);
 
436
g_assert(fp != NULL);
 
437
 
 
438
/* header */
 
439
fprintf(fp, "\n[ angles ]\n");
 
440
 
 
441
/* TODO - for all atoms with 2 + bonds -> check if 3 body term that matches */
 
442
for (list1=model->cores,k=1 ; list1 ; list1=g_slist_next(list1),k++)
 
443
  {
 
444
  core[1] = list1->data;
 
445
 
 
446
  nlist = connect_neighbours(core[1]);
 
447
 
 
448
  n = g_slist_length(nlist);
 
449
 
 
450
/* skip if too few, or suspiciously many */
 
451
  if (n < 2 || n >= MAX_CORES)
 
452
    continue;
 
453
 
 
454
  i = 0;
 
455
  for (list2 = nlist ; list2 ; list2=g_slist_next(list2))
 
456
    nc[i++] = list2->data;
 
457
 
 
458
/* iterate over unique neighbour pairs */
 
459
  for (i=0 ; i<n-1 ; i++)
 
460
    {
 
461
    core[0] = nc[i];
 
462
 
 
463
    for (j=i+1 ; j<n ; j++)
 
464
      {
 
465
      core[2] = nc[j];
 
466
 
 
467
      ff = ff_search(core, 3, model->ff_list);
 
468
 
 
469
      if (ff)
 
470
        {
 
471
/* TODO - write the entry */
 
472
/*
 
473
printf("FF found\n");
 
474
*/
 
475
 
 
476
/* get core numbering order */
 
477
        ni = g_slist_index(model->cores, core[0]) + 1;
 
478
        nj = g_slist_index(model->cores, core[2]) + 1;
 
479
 
 
480
        ff_gromacs = gromacs_type(ff);
 
481
 
 
482
/* lambda 0 */
 
483
/* print atoms and type */
 
484
        fprintf(fp, "%2d %2d %2d   %d  ", ni, k, nj, ff_gromacs->type);
 
485
 
 
486
/* NB: gromacs has the bond length 1st, followed by the potential parameters */
 
487
/* which is the reverse from internal/GULP format */
 
488
        if (ff_gromacs->bond_expected)
 
489
          fprintf(fp, "%f", ff_gromacs->bond_value);
 
490
        for (m=0 ; m<ff_gromacs->data_expected ; m++)
 
491
          fprintf(fp, " %f", ff_gromacs->data[m]);
 
492
 
 
493
/* lambda 1 */
 
494
        fprintf(fp, "    ");
 
495
        if (ff_gromacs->bond_expected)
 
496
          fprintf(fp, "%f", ff_gromacs->bond_value);
 
497
        for (m=0 ; m<ff_gromacs->data_expected ; m++)
 
498
          fprintf(fp, " %f", ff_gromacs->data[m]);
 
499
 
 
500
        fprintf(fp, "\n");
 
501
        }
 
502
      }
 
503
    }
 
504
  }
 
505
 
 
506
return(0);
 
507
}
 
508
 
 
509
/********************************/
 
510
/* write forcefield information */
 
511
/********************************/
 
512
gint write_gromacs_dihedrals(FILE *fp, struct model_pak *model)
 
513
{
 
514
gint i, j, k, l, n;
 
515
gdouble angle, da;
 
516
GSList *list, *list1, *list2, *end1, *end2;
 
517
GSList *d_list, *r_list, *i_list;
 
518
struct core_pak *c[4];
 
519
struct bond_pak *bond;
 
520
struct forcefield_pak *ff, *ff_gromacs;
 
521
 
 
522
g_assert(model != NULL);
 
523
 
 
524
d_list = ff_filter_list(FF_DIHEDRAL, model->ff_list);
 
525
r_list = ff_filter_list(FF_DIHEDRAL_RB, model->ff_list);
 
526
i_list = ff_filter_list(FF_IMPROPER, model->ff_list);
 
527
 
 
528
/* header */
 
529
fprintf(fp, "\n[ dihedrals ]\n");
 
530
 
 
531
/* generate torsions for all combinations about a bond */
 
532
for (list=model->bonds ; list ; list=g_slist_next(list))
 
533
  {
 
534
  bond = list->data;
 
535
 
 
536
  if (bond->type == BOND_HBOND)
 
537
    continue;
 
538
 
 
539
  c[1] = bond->atom1;
 
540
  c[2] = bond->atom2;
 
541
 
 
542
/* get lists of connected atoms at each end */
 
543
  end1 = connect_neighbours(c[1]);
 
544
  end1 = g_slist_remove(end1, c[2]);
 
545
  end2 = connect_neighbours(c[2]);
 
546
  end2 = g_slist_remove(end2, c[1]);
 
547
 
 
548
/* skip bonds with no neighbour atoms */
 
549
/*
 
550
  if (!end1 || !end2)
 
551
    goto gromacs_dihedral_loop_done;
 
552
*/
 
553
 
 
554
/* enumerate four body terms */
 
555
  for (list1=end1 ; list1 ; list1=g_slist_next(list1))
 
556
    {
 
557
    c[0] = list1->data;
 
558
 
 
559
    for (list2=end2 ; list2 ; list2=g_slist_next(list2))
 
560
      {
 
561
      c[3] = list2->data;
 
562
 
 
563
/* search for standard dihedral potential */
 
564
      ff = ff_search(c, 4, d_list);
 
565
 
 
566
/* search for RB dihedral */
 
567
      if (!ff)
 
568
        ff = ff_search(c, 4, r_list);
 
569
 
 
570
/* CURRENT */
 
571
/*
 
572
printf("[%s][%s][%s][%s]", c[0]->atom_label, c[1]->atom_label, c[2]->atom_label, c[3]->atom_label);
 
573
printf(" : dihedral = %f\n", measure_torsion(c));
 
574
*/
 
575
 
 
576
      if (ff)
 
577
        {
 
578
        ff_gromacs = gromacs_type(ff);
 
579
 
 
580
        i = g_slist_index(model->cores, c[0])+1;
 
581
        j = g_slist_index(model->cores, c[1])+1;
 
582
        k = g_slist_index(model->cores, c[2])+1;
 
583
        l = g_slist_index(model->cores, c[3])+1;
 
584
 
 
585
/* lambda 0 */
 
586
/* print atoms and type */
 
587
        fprintf(fp, "%2d %2d %2d %2d   %d  ", i, j, k, l, ff_gromacs->type);
 
588
/* NB: gromacs has the bond length 1st, followed by the potential parameters */
 
589
/* which is the reverse from internal/GULP format */
 
590
        if (ff_gromacs->bond_expected)
 
591
          fprintf(fp, "%f", ff_gromacs->bond_value);
 
592
        for (n=0 ; n<ff_gromacs->data_expected ; n++)
 
593
          fprintf(fp, " %f", ff_gromacs->data[n]);
 
594
 
 
595
/* lambda 1 */
 
596
        fprintf(fp, "    ");
 
597
        if (ff_gromacs->bond_expected)
 
598
          fprintf(fp, "%f", ff_gromacs->bond_value);
 
599
        for (n=0 ; n<ff_gromacs->data_expected ; n++)
 
600
          fprintf(fp, " %f", ff_gromacs->data[n]);
 
601
 
 
602
        fprintf(fp, "\n");
 
603
 
 
604
        g_free(ff_gromacs);
 
605
        }
 
606
      }
 
607
    }
 
608
 
 
609
  g_slist_free(end1);
 
610
  g_slist_free(end2);
 
611
  }
 
612
 
 
613
/* CURRENT - improper torsions (implemented as proper) */
 
614
/* FIXME - only valid for OPLS-AA */
 
615
/*
 
616
fprintf(fp, "\n[ impropers ]\n");
 
617
*/
 
618
fprintf(fp, "; Improper terms, implemented as proper torsions.\n");
 
619
 
 
620
/* for each atom with 3 bonds (planar) */
 
621
for (list=model->cores ; list ; list=g_slist_next(list))
 
622
  {
 
623
  c[0] = list->data;
 
624
 
 
625
  list1 = connect_neighbours(c[0]); 
 
626
 
 
627
/* CHECK - only an atom with 3 bonds can be a candidate */
 
628
  if (g_slist_length(list1) == 3)
 
629
    {
 
630
    c[1] = list1->data;
 
631
    list2 = g_slist_next(list1);
 
632
    c[2] = list2->data;
 
633
    list2 = g_slist_next(list2);
 
634
    c[3] = list2->data;
 
635
 
 
636
 
 
637
/* TODO - want c[3] to be defined so that the c[1] - c[2] axis is */
 
638
/* orthogonal to the c[0] - c[3] direction */
 
639
 
 
640
/* TODO - choose c[1] & c[2] to be same type (if possible) */
 
641
 
 
642
/* improper torsion - only allow if angle is ~180 */
 
643
    angle = measure_torsion(c);
 
644
/*
 
645
    da = fabs(180.0 - fabs(angle));
 
646
*/
 
647
    da = fabs(angle);
 
648
 
 
649
/* 10 degree tolerance */
 
650
    ff = NULL;
 
651
    if (da < 10.0)
 
652
      ff = ff_search(c, 4, i_list);
 
653
 
 
654
    if (ff)
 
655
      {
 
656
      ff_gromacs = gromacs_type(ff);
 
657
 
 
658
      i = g_slist_index(model->cores, c[0])+1;
 
659
      j = g_slist_index(model->cores, c[1])+1;
 
660
      k = g_slist_index(model->cores, c[2])+1;
 
661
      l = g_slist_index(model->cores, c[3])+1;
 
662
 
 
663
/* lambda 0 */
 
664
      fprintf(fp, "%2d %2d %2d %2d   1 ", i, j, k, l);
 
665
      fprintf(fp, " %f", ff_gromacs->bond_value);
 
666
      for (n=0 ; n<ff_gromacs->data_current ; n++)
 
667
        fprintf(fp, " %f", ff_gromacs->data[n]);
 
668
 
 
669
/* lambda 1 */
 
670
      fprintf(fp, "     %f", ff_gromacs->bond_value);
 
671
      for (n=0 ; n<ff_gromacs->data_current ; n++)
 
672
        fprintf(fp, " %f", ff_gromacs->data[n]);
 
673
 
 
674
      fprintf(fp, "\n");
 
675
 
 
676
      g_free(ff_gromacs);
 
677
      }
 
678
    }
 
679
  g_slist_free(list1);
 
680
  }
 
681
 
 
682
/* cleanup */
 
683
g_slist_free(d_list);
 
684
g_slist_free(r_list);
 
685
g_slist_free(i_list);
 
686
 
 
687
return(0);
 
688
}
 
689
 
 
690
/********************************/
 
691
/* write forcefield information */
 
692
/********************************/
 
693
gint write_gromacs_nonbond(FILE *fp, struct model_pak *model)
 
694
{
 
695
gint i;
 
696
GSList *list;
 
697
struct forcefield_pak *ff, *ff_gromacs;
 
698
 
 
699
g_assert(model != NULL);
 
700
 
 
701
/* header */
 
702
fprintf(fp, "\n[ nonbond_params ]\n");
 
703
 
 
704
for (list=model->ff_list ; list ; list=g_slist_next(list))
 
705
  {
 
706
  ff = list->data;
 
707
 
 
708
/* skip anything that doesnt have only 2 atom identifiers */
 
709
  if (ff->atoms_expected != 2)
 
710
    continue;
 
711
  if (ff->atoms_current != 2)
 
712
    continue;
 
713
 
 
714
/* GROMACS only supports LJ/BUCK */
 
715
  if (ff->type == FF_LENNARD || ff->type == FF_BUCKINGHAM)
 
716
    {
 
717
/* get appropriate type */
 
718
    ff_gromacs = gromacs_type(ff);
 
719
 
 
720
/* print LJ/Buck */
 
721
    fprintf(fp, "%s %s %d", ff->atom[0], ff->atom[1], ff_gromacs->type);
 
722
    for (i=0 ; i<ff_gromacs->data_expected ; i++)
 
723
      fprintf(fp, " %f", ff_gromacs->data[i]);
 
724
    fprintf(fp, "\n");
 
725
 
 
726
    }
 
727
  }
 
728
 
 
729
return(0);
 
730
}
 
731
 
 
732
/***********************************************************/
 
733
/* topology ie molecule, connectivity, force fields (.top) */
 
734
/***********************************************************/
 
735
gint write_gromacs_top(gchar *filename, struct model_pak *model)
 
736
{
 
737
FILE *fp;
 
738
 
 
739
fp = fopen(filename, "wt");
 
740
 
 
741
/* for each molecule write topology entry */
 
742
/* TODO - sort into unique molecules (atoms must also be ordered the same) */
 
743
 
 
744
write_gromacs_header(fp, model);
 
745
 
 
746
write_gromacs_atoms(fp, model);
 
747
 
 
748
write_gromacs_bonds(fp, model);
 
749
 
 
750
write_gromacs_angles(fp, model);
 
751
 
 
752
write_gromacs_dihedrals(fp, model);
 
753
 
 
754
write_gromacs_nonbond(fp, model);
 
755
 
 
756
write_gromacs_footer(fp, model);
 
757
 
 
758
fclose(fp);
 
759
 
 
760
return(0);
 
761
}
 
762
 
 
763
/****************/
 
764
/* file writing */
 
765
/****************/
 
766
gint write_gromacs(gchar *filename, struct model_pak *model)
 
767
{
 
768
gchar *temp;
 
769
 
 
770
temp = parse_extension_set(filename, "top");
 
771
 
 
772
write_gromacs_gro(filename, model);
 
773
write_gromacs_top(temp, model);
 
774
 
 
775
g_free(temp);
 
776
 
 
777
return(0);
 
778
}
 
779
 
 
780
/********************************/
 
781
/* read in a GROMACS FF library */
 
782
/********************************/
 
783
#define DEBUG_GROMACS_READ_FF 0
 
784
GSList *gromacs_read_ff(const gchar *filename)
 
785
{
 
786
gint state, type, num_tokens, ub, ua, ud, i;
 
787
gchar *line, **buff, **abuff, **dbuff;
 
788
GSList *list=NULL;
 
789
struct forcefield_pak *ff;
 
790
FILE *fp;
 
791
 
 
792
fp = fopen(filename, "rt");
 
793
if (!fp)
 
794
  return(NULL);
 
795
 
 
796
/* unknown bond, angle, dihedral types */
 
797
ub = ua = ud = 0;
 
798
 
 
799
/* stop auto skip of lines starting with # ... ie #define */
 
800
file_skip_comment(FALSE);
 
801
 
 
802
line = file_read_line(fp);
 
803
state = -1;
 
804
while (line)
 
805
  {
 
806
/* TODO - get rid of leading whitespace, so funny comment lines are properly ignored */
 
807
 
 
808
  ff = NULL;
 
809
 
 
810
/* directive processing */
 
811
  switch(line[0])
 
812
    {
 
813
    case '[':
 
814
      state = -1;
 
815
      if (g_strrstr(line, "bond"))
 
816
        state = BOND;
 
817
      if (g_strrstr(line, "angle"))
 
818
        state = ANGLE;
 
819
      if (g_strrstr(line, "dihedral"))
 
820
        state = DIHEDRAL;
 
821
      break;
 
822
 
 
823
/* process #defines */
 
824
    case '#':
 
825
      if (g_ascii_strncasecmp("#define", line, 7) == 0) 
 
826
        {
 
827
        buff = get_tokens(line, 3);
 
828
        abuff = g_strsplit(*(buff+1), "_", -1);
 
829
 
 
830
#if GTK_MAJOR_VERSION >= 2 && GTK_MINOR_VERSION >= 6
 
831
        num_tokens = g_strv_length(abuff);
 
832
#else
 
833
        num_tokens = 0;
 
834
#endif
 
835
 
 
836
        if (abuff)
 
837
          {
 
838
          if (g_ascii_strncasecmp(*(abuff), "improper", 8) == 0)
 
839
            {
 
840
            if (num_tokens > 4)
 
841
              {
 
842
              ff = ff_type_new(FF_IMPROPER);
 
843
              g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(abuff+num_tokens-4));
 
844
              g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(abuff+num_tokens-3));
 
845
              g_snprintf(ff->atom[2], FF_MAX_SYMBOL, "%s", *(abuff+num_tokens-2));
 
846
              g_snprintf(ff->atom[3], FF_MAX_SYMBOL, "%s", *(abuff+num_tokens-1));
 
847
/* TODO - Y, Z -> X ... since X is wildcard type */
 
848
              for (i=4 ; i-- ; )
 
849
                if (ff->atom[i][0] == 'Y' || ff->atom[i][0] == 'Z')
 
850
                  ff->atom[i][0] = 'X';
 
851
 
 
852
              dbuff = get_tokens(line, 6);
 
853
              ff->bond_value = str_to_float(*(dbuff+2));
 
854
              ff->bond_units = FF_DEG;
 
855
              ff->data[0] = str_to_float(*(dbuff+3));
 
856
              ff->data[1] = str_to_float(*(dbuff+4));
 
857
              g_strfreev(dbuff);
 
858
 
 
859
              ff->atoms_current = ff->atoms_expected;
 
860
              ff->data_current = ff->data_expected = 2;
 
861
 
 
862
              list = g_slist_prepend(list, ff);
 
863
              ff = NULL;
 
864
              }
 
865
            }
 
866
          g_strfreev(abuff);
 
867
          }
 
868
        g_strfreev(buff);
 
869
        }
 
870
      state = -1;
 
871
      break;
 
872
 
 
873
/* comment */
 
874
    case ';':
 
875
      break;
 
876
 
 
877
    default:
 
878
/* normal processing */
 
879
      buff = tokenize(line, &num_tokens);
 
880
      switch (state)
 
881
        {
 
882
        case BOND:
 
883
          if (num_tokens > 4)
 
884
            {
 
885
            type = g_ascii_strtod(*(buff+2), NULL);
 
886
            switch (type)
 
887
              {
 
888
              case 1:
 
889
                ff = ff_type_new(FF_HARMONIC);
 
890
                g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(buff+0));
 
891
                g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(buff+1));
 
892
/* NB: GROMACS uses nm rather than angs */
 
893
                ff->bond_value = 10.0 * str_to_float(*(buff+3));
 
894
                ff->bond_units = FF_ANG;
 
895
/* nm -> ang correction */
 
896
                ff->data[0] = 0.01 * str_to_float(*(buff+4));
 
897
                ff->data_units = FF_KJ;
 
898
 
 
899
                ff->atoms_current = ff->atoms_expected;
 
900
                ff->data_current = ff->data_expected;
 
901
                break;
 
902
 
 
903
              case 3:
 
904
g_assert(num_tokens > 5);
 
905
                ff = ff_type_new(FF_MORSE);
 
906
                g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(buff+0));
 
907
                g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(buff+1));
 
908
/* nm -> ang correction */
 
909
                ff->bond_value = 10.0*str_to_float(*(buff+3));
 
910
                ff->bond_units = FF_ANG;
 
911
                ff->data[0] = str_to_float(*(buff+4));
 
912
                ff->data[1] = 0.1 * str_to_float(*(buff+5));
 
913
                ff->data_units = FF_KJ;
 
914
 
 
915
                ff->atoms_current = ff->atoms_expected;
 
916
                ff->data_current = ff->data_expected;
 
917
                break;
 
918
 
 
919
              default:
 
920
                ub++;
 
921
              }
 
922
            }
 
923
          break;
 
924
 
 
925
        case ANGLE:
 
926
          if (num_tokens > 5)
 
927
            {
 
928
            type = g_ascii_strtod(*(buff+3), NULL);
 
929
            switch (type)
 
930
              {
 
931
              case 1:
 
932
                ff = ff_type_new(FF_3B_HARMONIC);
 
933
                g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(buff+0));
 
934
                g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(buff+1));
 
935
                g_snprintf(ff->atom[2], FF_MAX_SYMBOL, "%s", *(buff+2));
 
936
                ff->bond_value = str_to_float(*(buff+4));
 
937
                ff->bond_units = FF_DEG;
 
938
                ff->data[0] = str_to_float(*(buff+5));
 
939
/* strictly - kJ/mol rad-2 */
 
940
                ff->data_units = FF_KJ;
 
941
                ff->atoms_current = ff->atoms_expected;
 
942
                ff->data_current = ff->data_expected;
 
943
                break;
 
944
 
 
945
              default:
 
946
                ua++;
 
947
              }
 
948
            }
 
949
          break;
 
950
 
 
951
        case DIHEDRAL:
 
952
          if (num_tokens > 4)
 
953
            {
 
954
            type = g_ascii_strtod(*(buff+4), NULL);
 
955
            switch (type)
 
956
              {
 
957
              case 1:
 
958
                if (num_tokens < 7)
 
959
                  {
 
960
                  ud++;
 
961
                  break;
 
962
                  }
 
963
                ff = ff_type_new(FF_DIHEDRAL);
 
964
                g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(buff+0));
 
965
                g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(buff+1));
 
966
                g_snprintf(ff->atom[2], FF_MAX_SYMBOL, "%s", *(buff+2));
 
967
                g_snprintf(ff->atom[3], FF_MAX_SYMBOL, "%s", *(buff+3));
 
968
                ff->bond_value = str_to_float(*(buff+5));
 
969
                ff->bond_units = FF_DEG;
 
970
                ff->data[0] = str_to_float(*(buff+6));
 
971
                ff->data_units = FF_KJ;
 
972
                ff->atoms_current = ff->atoms_expected;
 
973
                ff->data_current = ff->data_expected;
 
974
                break;
 
975
 
 
976
              case 3:
 
977
                if (num_tokens < 11)
 
978
                  {
 
979
                  ud++;
 
980
                  break;
 
981
                  }
 
982
                ff = ff_type_new(FF_DIHEDRAL_RB);
 
983
                g_snprintf(ff->atom[0], FF_MAX_SYMBOL, "%s", *(buff+0));
 
984
                g_snprintf(ff->atom[1], FF_MAX_SYMBOL, "%s", *(buff+1));
 
985
                g_snprintf(ff->atom[2], FF_MAX_SYMBOL, "%s", *(buff+2));
 
986
                g_snprintf(ff->atom[3], FF_MAX_SYMBOL, "%s", *(buff+3));
 
987
                ff->data[0] = str_to_float(*(buff+5));
 
988
                ff->data[1] = str_to_float(*(buff+6));
 
989
                ff->data[2] = str_to_float(*(buff+7));
 
990
                ff->data[3] = str_to_float(*(buff+8));
 
991
                ff->data[4] = str_to_float(*(buff+9));
 
992
                ff->data[5] = str_to_float(*(buff+10));
 
993
                ff->data_units = FF_KJ;
 
994
 
 
995
                ff->atoms_current = ff->atoms_expected;
 
996
                ff->data_current = ff->data_expected;
 
997
                break;
 
998
 
 
999
              default:
 
1000
                ud++;
 
1001
              }
 
1002
            }
 
1003
          break;
 
1004
        }
 
1005
 
 
1006
/* add the object to the return list */
 
1007
      if (ff)
 
1008
        list = g_slist_prepend(list, ff);
 
1009
 
 
1010
      g_strfreev(buff);
 
1011
      break;
 
1012
    }
 
1013
 
 
1014
  g_free(line);
 
1015
  line = file_read_line(fp);
 
1016
  }
 
1017
 
 
1018
g_free(line);
 
1019
 
 
1020
#if DEBUG_GROMACS_READ_FF
 
1021
printf("processed: %d items.\n", g_slist_length(list));
 
1022
printf("ignored: [bonds = %d] [angles = %d] [dihedrals = %d]\n", ub, ua, ud);
 
1023
/*
 
1024
ff_dump_type(FF_IMPROPER, list);
 
1025
*/
 
1026
#endif
 
1027
 
 
1028
file_skip_comment(TRUE);
 
1029
 
 
1030
return(list);
 
1031
}