?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix compile error for IBM VMX on ppc64
[gromacs/AngularHB.git] / src / contrib / mkice.c
bloba60b847cdd9b881be02bcde7221ec663c71cf627
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.3.99_development_20071104
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-2006, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Groningen Machine for Chemical Simulation
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <stdio.h>
40 #include <math.h>
41 #include "typedefs.h"
42 #include "gromacs/commandline/pargs.h"
43 #include "copyrite.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/fileio/pdbio.h"
46 #include "macros.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/pbcutil/pbc.h"
50 #include "gromacs/math/units.h"
51 #include "names.h"
52 #include "txtdump.h"
53 #include "gromacs/fileio/trnio.h"
54 #include "gromacs/topology/symtab.h"
55 #include "gromacs/fileio/strdb.h"
56 #include "gromacs/fileio/confio.h"
57
58 #define TET 109.47
59 #define DCONS 0.117265878
60
61 typedef struct {
62 int n,aa[4];
63 } t_bbb;
64
65 static char *watname[] = { "OW ", "HW1", "HW2", "DW", "SW" };
66 static char *diamname[] = { "C1", "H2", "H3", "H4", "H5", "H2", "H3", "H4", "H5" };
67 static real qspc[] = { -0.82, 0.41, 0.41 };
68 static real qyaw[] = { 1.24588, 0.62134, 0.62134, 0.0, -2.48856 };
69 static real spc_lj[3][6] = {
70 { 2.6171e-3, 2.6331e-6, 0, 0, 0, 0 },
71 { 0, 0, 0, 0, 0, 0 },
72 { 0, 0, 0, 0, 0, 0 }
73 };
74 #define CHH 2e-9
75 #define CHS 2e-9
76 #define COS 2e-6
77 static real yaw_lj[5][10] = {
78 { 0, 0, 0, 0, 0, 0, 0, 0, 0, COS },
79 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
80 { 0, 0, 0, CHH, 0, CHH, 0, 0, 0, CHS },
81 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
82 { 0, COS, 0, CHS, 0, CHS, 0, 0, 2.6e-3, 0 }
83 };
84
85 void unitcell(rvec x[],rvec box,gmx_bool bYaw,real odist,real hdist)
86 {
87 #define cx 0.81649658
88 #define cy 0.47140452
89 #define cy2 0.94280904
90 #define cz 0.33333333
91
92 rvec xx[24] = {
93 { 0, 0, 0 }, /* O1 */
94 { 0, 0, 1 }, /* H relative to Oxygen */
95 { cx, cy, -cz },
96 { cx, cy, -cz }, /* O2 */
97 { 0, 0, -1 }, /* H relative to Oxygen */
98 { cx,-cy, +cz },
99 { cx, cy+cy2, 0 }, /* O3 */
100 { -cx, cy, -cz }, /* H relative to Oxygen */
101 { 0, -cy2, -cz },
102 { 0, 2*cy+cy2, -cz }, /* O4 */
103 {-cx,-cy, +cz }, /* H relative to Oxygen */
104 { 0 , cy2, +cz },
105 { 0, 0, 1 }, /* O5 */
106 {-cx, cy, +cz }, /* H relative to Oxygen */
107 { 0 , -cy2, +cz },
108 { cx, cy, 1+cz }, /* O6 */
109 { -cx, -cy, -cz }, /* H relative to Oxygen */
110 { 0, cy2, -cz },
111 { cx, cy+cy2, 1 }, /* O7 */
112 { 0, 0, -1 }, /* H relative to Oxygen */
113 { cx, cy, +cz },
114 { 0, 2*cy+cy2,1+cz }, /* O8 */
115 { 0, 0, 1 }, /* H relative to Oxygen */
116 { cx, -cy, -cz }
117 };
118 int i,iin,iout,j,m;
119 rvec tmp,t2,dip;
120
121 clear_rvec(dip);
122 for(i=0; (i<8); i++) {
123 iin = 3*i;
124 if (bYaw)
125 iout = 5*i;
126 else
127 iout = iin;
128 svmul(odist,xx[iin],x[iout]);
129 svmul(-0.82,x[iout],t2);
130 rvec_inc(dip,t2);
131 for(j=1; (j<=2); j++) {
132 svmul(hdist,xx[iin+j],tmp);
133 rvec_add(x[iout],tmp,x[iout+j]);
134 svmul(0.41,x[iout+j],t2);
135 rvec_inc(dip,t2);
136 }
137 if (bYaw)
138 for(m=0; (m<DIM); m++)
139 x[iout+3][m] = x[iout+4][m] =
140 (1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
141 }
142
143 box[XX] = 2*cx;
144 box[YY] = 2*(cy2+cy);
145 box[ZZ] = 2*(1+cz);
146 for(i=0; (i<DIM); i++)
147 box[i] *= odist;
148
149 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
150 printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
151 }
152
153 void random_h_coords(int natmol,int nmol,rvec x[],rvec box,
154 gmx_bool bYaw,real odist,real hdist)
155 {
156 #define cx 0.81649658
157 #define cy 0.47140452
158 #define cy2 0.94280904
159 #define cz 0.33333333
160
161 rvec xx[24] = {
162 { 0, 0, 0 }, /* O1 */
163 { 0, 0, 1 }, /* H relative to Oxygen */
164 { cx, cy, -cz },
165 { cx, cy, -cz }, /* O2 */
166 { 0, 0, -1 }, /* H relative to Oxygen */
167 { cx,-cy, +cz },
168 { cx, cy+cy2, 0 }, /* O3 */
169 { -cx, cy, -cz }, /* H relative to Oxygen */
170 { 0, -cy2, -cz },
171 { 0, 2*cy+cy2, -cz }, /* O4 */
172 {-cx,-cy, +cz }, /* H relative to Oxygen */
173 { 0 , cy2, +cz },
174 { 0, 0, 1 }, /* O5 */
175 {-cx, cy, +cz }, /* H relative to Oxygen */
176 { 0 , -cy2, +cz },
177 { cx, cy, 1+cz }, /* O6 */
178 { -cx, -cy, -cz }, /* H relative to Oxygen */
179 { 0, cy2, -cz },
180 { cx, cy+cy2, 1 }, /* O7 */
181 { 0, 0, -1 }, /* H relative to Oxygen */
182 { cx, cy, +cz },
183 { 0, 2*cy+cy2,1+cz }, /* O8 */
184 { 0, 0, 1 }, /* H relative to Oxygen */
185 { cx, -cy, -cz }
186 };
187 int i,iin,iout,j,m;
188 rvec tmp,t2,dip;
189
190 clear_rvec(dip);
191 for(i=0; (i<nmol); i++) {
192 iin = natmol*i;
193 iout = iin;
194 svmul(odist,x[iin],x[iout]);
195 svmul(-0.82,x[iout],t2);
196 rvec_inc(dip,t2);
197 for(j=1; (j<=2); j++) {
198 svmul(hdist,xx[3*(i % 8)+j],tmp);
199 rvec_add(x[iout],tmp,x[iout+j]);
200 svmul(0.41,x[iout+j],t2);
201 rvec_inc(dip,t2);
202 }
203 }
204
205 box[XX] = 2*cx;
206 box[YY] = 2*(cy2+cy);
207 box[ZZ] = 2*(1+cz);
208 for(i=0; (i<DIM); i++)
209 box[i] *= odist;
210
211 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
212 printf("Dipole: %10.5f %10.5f %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
213 }
214
215 void unitcell_d(rvec x[],rvec box,real odist)
216 {
217 rvec cc[8] = {
218 { 0, 0, 0 }, /* C1 */
219 { cx, cy, -cz }, /* C2 */
220 { cx, cy+cy2, 0 }, /* C3 */
221 { 0, 2*cy+cy2, -cz }, /* C4 */
222 { 0, 0, 1 }, /* C5 */
223 { cx, cy, 1+cz }, /* C6 */
224 { cx, cy+cy2, 1 }, /* C7 */
225 { 0, 2*cy+cy2,1+cz }, /* C8 */
226 };
227 rvec hh[4] = {
228 { 0, 0, 1 }, /* H relative to C */
229 { cx, cy, -cz },
230 { cx, -cy, -cz },
231 {-cy2, 0, -cz }
232 };
233 int i,iin,iout,j,m;
234 rvec tmp,t2,dip;
235
236 clear_rvec(dip);
237 for(i=0; (i<8); i++) {
238 iin = i;
239 iout = i;
240 svmul(odist,cc[iin],x[iout]);
241 }
242 box[XX] = 2*cx;
243 box[YY] = 2*(cy2+cy);
244 box[ZZ] = 2*(1+cz);
245 for(i=0; (i<DIM); i++)
246 box[i] *= odist;
247
248 printf("Unitcell: %10.5f %10.5f %10.5f\n",box[XX],box[YY],box[ZZ]);
249 }
250
251 static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
252 gmx_bool bPBC,matrix box)
253 {
254 real od2 = odist*odist+1e-5;
255 t_pbc pbc;
256 t_bbb *bbb;
257 int i,j;
258 rvec dx;
259
260 if (bPBC)
261 set_pbc(&pbc,box);
262 snew(bbb,natoms);
263 for(i=0; (i<natoms); i++) {
264 for(j=i+1; (j<natoms); j++) {
265 if (bPBC)
266 pbc_dx(&pbc,x[i],x[j],dx);
267 else
268 rvec_sub(x[i],x[j],dx);
269 if (iprod(dx,dx) <= od2) {
270 bbb[i].aa[bbb[i].n++] = j;
271 bbb[j].aa[bbb[j].n++] = i;
272 }
273 }
274 }
275 if (debug)
276 #define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
277 for(i=0; (i<natoms); i++)
278 fprintf(debug,"bonds from %3d: %d %d %d %d\n",
279 i,PRB(1),PRB(2),PRB(3),PRB(4));
280 #undef PRB
281 return bbb;
282 }
283
284 static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
285 gmx_bool bPBC,matrix box)
286 {
287
288 int i,ib,j,k,l,m,nrm=0;
289 t_bbb *bbb;
290 gmx_bool *bRemove;
291 rvec dx;
292
293 do {
294 nrm = 0;
295 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
296
297 for(i=0; (i<a->nr); i++) {
298 if (bbb[i].n < 2) {
299 for(k=0; (k<bbb[i].n); k++) {
300 ib = bbb[i].aa[k];
301 for(j=0; (j<bbb[ib].n); j++)
302 if (bbb[ib].aa[j] == i)
303 break;
304 if (j == bbb[ib].n)
305 gmx_fatal(FARGS,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
306 for( ; (j<bbb[ib].n-1); j++)
307 bbb[ib].aa[j] = bbb[ib].aa[j+1];
308 bbb[ib].n--;
309 nrm++;
310 }
311 bbb[i].n = 0;
312 }
313 }
314
315 for(i=j=0; (i<a->nr); i++) {
316 if (bbb[i].n >= 2) {
317 copy_rvec(x[i],x[j]);
318 j++;
319 }
320 }
321 fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
322 a->nr-j,a->nr);
323 a->nr = j;
324 sfree(bbb);
325 } while (nrm > 0);
326 /* Rename atoms */
327 bbb = mk_bonds(a->nr,x,odist,bPBC,box);
328 for(i=0; (i<a->nr); i++) {
329 switch (bbb[i].n) {
330 case 4:
331 a->atomname[i] = put_symtab(symtab,"C");
332 break;
333 case 3:
334 a->atomname[i] = put_symtab(symtab,"CH1");
335 break;
336 case 2:
337 a->atomname[i] = put_symtab(symtab,"CH2");
338 break;
339 default:
340 gmx_fatal(FARGS,"This atom (%d) has %d bonds only",i,bbb[i].n);
341 }
342 }
343 sfree(bbb);
344 }
345
346 static real calc_ener(real c6,real c12,rvec dx,tensor vir)
347 {
348 real r,e,f;
349 int m,n;
350
351 r = norm(dx);
352 e = c12*pow(r,-12.0) - c6*pow(r,-6.0);
353 f = 12*c12*pow(r,-14.0) - 6*c6*pow(r,-8.0);
354 for(m=0; (m<DIM); m++)
355 for(n=0; (n<DIM); n++)
356 vir[m][n] -= 0.5*f*dx[m]*dx[n];
357
358 return e;
359 }
360
361 static int read_rel_coords(char *fn,rvec **xx,int natmol)
362 {
363 int i,nline;
364 double x,y,z;
365 char **strings=NULL;
366
367 nline = get_file(fn,&strings);
368 printf("Read %d lines from %s\n",nline,fn);
369 snew(*xx,nline*natmol);
370 for(i=0; (i<nline); i++) {
371 if (sscanf(strings[i],"%lf%lf%lf",&x,&y,&z) != 3)
372 gmx_fatal(FARGS,"Not enough arguments on line %d of file %s (should be 3)",i,fn);
373 (*xx)[natmol*i][XX] = x;
374 (*xx)[natmol*i][YY] = y;
375 (*xx)[natmol*i][ZZ] = z;
376 }
377 return natmol*nline;
378 }
379
380 void virial(FILE *fp,gmx_bool bFull,int nmol,rvec x[],matrix box,real rcut,
381 gmx_bool bYaw,real q[],gmx_bool bLJ)
382 {
383 int i,j,im,jm,natmol,ik,jk,m,ninter;
384 rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
385 real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
386 t_pbc pbc;
387
388 set_pbc(&pbc,box);
389 fprintf(fp,"%3s - %3s: %6s %6s %6s %6s %8s %8s %8s\n",
390 "ai","aj","dx","dy","dz","|d|","virx","viry","virz");
391 clear_rvec(ftot);
392 clear_rvec(vir);
393 ninter = 0;
394 vctot = 0;
395 v12tot = 0;
396 v6tot = 0;
397 natmol = bYaw ? 5 : 3;
398 snew(force,nmol*natmol);
399
400 for(i=0; (i<nmol); i++) {
401 im = natmol*i;
402 /* Center of geometry */
403 clear_rvec(xcmi);
404 for(ik=0; (ik<natmol); ik++)
405 rvec_inc(xcmi,x[im+ik]);
406 for(m=0; (m<DIM); m++)
407 xcmi[m] /= natmol;
408
409 for(j=i+1; (j<nmol); j++) {
410 jm = natmol*j;
411 /* Center of geometry */
412 clear_rvec(xcmj);
413 for(jk=0; (jk<natmol); jk++)
414 rvec_inc(xcmj,x[jm+jk]);
415 for(m=0; (m<DIM); m++)
416 xcmj[m] /= natmol;
417
418 /* First check COM-COM distance */
419 pbc_dx(&pbc,xcmi,xcmj,dx);
420 if (norm(dx) < rcut) {
421 ninter++;
422 /* Neirest neighbour molecules! */
423 clear_rvec(dvir);
424 for(ik=0; (ik<natmol); ik++) {
425 for(jk=0; (jk<natmol); jk++) {
426 pbc_dx(&pbc,x[im+ik],x[jm+jk],dx);
427 dx2 = iprod(dx,dx);
428 dx1 = sqrt(dx2);
429 vcoul = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
430 vctot += vcoul;
431
432 if (bLJ) {
433 if (bYaw) {
434 c6 = yaw_lj[ik][2*jk];
435 c12 = yaw_lj[ik][2*jk+1];
436 }
437 else {
438 c6 = spc_lj[ik][2*jk];
439 c12 = spc_lj[ik][2*jk+1];
440 }
441 dx6 = dx2*dx2*dx2;
442 v6 = c6/dx6;
443 v12 = c12/(dx6*dx6);
444 v6tot -= v6;
445 v12tot+= v12;
446 fscal = (vcoul+12*v12-6*v6)/dx2;
447 }
448 else
449 fscal = vcoul/dx2;
450 for(m=0; (m<DIM); m++) {
451 f[m] = dx[m]*fscal;
452 dvir[m] -= 0.5*dx[m]*f[m];
453 }
454 rvec_inc(force[ik+im],f);
455 rvec_dec(force[jk+jm],f);
456 /*if (bFull)
457 fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
458 " %8.3f %8.3f %8.3f\n",
459 watname[ik],im+ik,watname[jk],jm+jk,
460 dx[XX],dx[YY],dx[ZZ],norm(dx),
461 dvir[XX],dvir[YY],dvir[ZZ]);*/
462 }
463 }
464 if (bFull)
465 fprintf(fp,"%3s%4d-%3s%4d: "
466 " %8.3f %8.3f %8.3f\n",
467 "SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
468 rvec_inc(vir,dvir);
469 }
470 }
471 }
472 fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
473 ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
474 fprintf(fp,"Vcoul: %10.4e V12: %10.4e V6: %10.4e Vtot: %10.4e (kJ/mol)\n",
475 vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
476 pr_rvec(fp,0,"vir ",vir,DIM,TRUE);
477
478 for(m=0; (m<DIM); m++)
479 pres[m] = -2*PRESFAC/(det(box))*vir[m];
480 pr_rvec(fp,0,"pres",pres,DIM,TRUE);
481 pr_rvecs(fp,0,"force",force,natmol*nmol);
482 sfree(force);
483 }
484
485
486
487 int main(int argc,char *argv[])
488 {
489 static char *desc[] = {
490 "[TT]mkice[tt] generates an ice crystal in the Ih crystal form which is the",
491 "most stable form. The rectangular unitcell contains eight molecules",
492 "and all oxygens are tetrahedrally coordinated.[PAR]",
493 "If an input file is given it is interpreted as a series of oxygen",
494 "coordinates the distance between which can be scaled by the odist flag.",
495 "The program then adds hydrogens to the oxygens in random orientation",
496 "but with proper OH distances and HOH angle. This feature allows to",
497 "build water clusters based on oxygen coordinates only."
498 };
499 static int nx=1,ny=1,nz=1;
500 static gmx_bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE;
501 static gmx_bool bOrdered=TRUE,bDiamond=FALSE,bPBC=TRUE;
502 static real rcut=0.3,odist=0.274,hdist=0.09572;
503 t_pargs pa[] = {
504 { "-nx", FALSE, etINT, {&nx}, "nx" },
505 { "-ny", FALSE, etINT, {&ny}, "ny" },
506 { "-nz", FALSE, etINT, {&nz}, "nz" },
507 { "-yaw", FALSE, etBOOL, {&bYaw},
508 "Generate virtual sites and shell positions" },
509 { "-lj", FALSE, etBOOL, {&bLJ},
510 "Use LJ as well as coulomb for virial calculation" },
511 { "-rcut", FALSE,etREAL, {&rcut},"Cut-off for virial calculations" },
512 { "-full", FALSE,etBOOL, {&bFull},"Full virial output" },
513 { "-odist", FALSE, etREAL, {&odist}, "Distance (nm) between oxygens" },
514 { "-hdist", FALSE, etREAL, {&hdist}, "Bondlength (nm) for OH bond" },
515 { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
516 { "-pbc", FALSE, etBOOL, {&bPBC}, "Make a periodic diamond" },
517 { "-order", FALSE,etBOOL, {&bOrdered}, "Make a proton-ordered ice lattice" },
518 { "-series",FALSE, etBOOL, {&bSeries},
519 "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
520 };
521 t_filenm fnm[] = {
522 { efSTO, "-p", "ice", ffWRITE },
523 { efSTO, "-c", NULL, ffOPTRD },
524 { efDAT, "-f", NULL, ffOPTRD },
525 { efTRN, "-o", "ice", ffOPTWR }
526 };
527 #define NFILE asize(fnm)
528
529 FILE *fp;
530 char *fn,quote[256];
531 int i,j,k,n,nmax,m,natom,natmol;
532 t_atoms *pdba;
533 t_atoms atoms;
534 t_symtab symtab;
535 rvec box,tmp,*xx;
536 matrix boxje;
537
538 CopyRight(stdout,argv[0]);
539 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),
540 desc,0,NULL);
541 if (debug) {
542 fprintf(debug,"nx = %3d, ny = %3d, nz = %3d\n",nx,ny,nz);
543 fprintf(debug,"YAW = %3s, LJ = %3s, rcut = %g\n",yesno_names[bYaw],
544 yesno_names[bLJ],rcut);
545 }
546
547 if (bYaw)
548 natmol = 5;
549 else if (bDiamond)
550 natmol = 1;
551 else
552 natmol = 3;
553
554 if (opt2bSet("-f",NFILE,fnm)) {
555 natom = read_rel_coords(opt2fn("-f",NFILE,fnm),&xx,natmol);
556 nmax = natom;
557 }
558 else {
559 natom = natmol*8;
560 nmax = natom*nx*ny*nz;
561 snew(xx,nmax);
562 }
563 snew(pdba,1);
564 init_t_atoms(pdba,nmax,TRUE);
565 pdba->nr = nmax;
566 open_symtab(&symtab);
567 for(n=0; (n<nmax); n++) {
568 pdba->pdbinfo[n].type = epdbATOM;
569 pdba->pdbinfo[n].atomnr = 1+n;
570 pdba->atom[n].resnr = 1+(n/natmol);
571 pdba->atomname[n] = put_symtab(&symtab,
572 bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
573 if (bDiamond)
574 pdba->resname[n] = put_symtab(&symtab,"DIA");
575 else
576 pdba->resname[n] = put_symtab(&symtab,"SOL");
577 sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
578 pdba->atom[n].chain = ' ';
579 }
580
581 /* Generate the unit cell */
582 if (bDiamond)
583 unitcell_d(xx,box,odist);
584 else if (opt2bSet("-f",NFILE,fnm)) {
585 random_h_coords(natmol,natom/natmol,xx,box,bYaw,odist,hdist);
586 }
587 else
588 unitcell(xx,box,bYaw,odist,hdist);
589 if (debug) {
590 clear_mat(boxje);
591 boxje[XX][XX] = box[XX];
592 boxje[YY][YY] = box[YY];
593 boxje[ZZ][ZZ] = box[ZZ];
594 }
595 n=0;
596 for(i=0; (i<nx); i++) {
597 tmp[XX] = i*box[XX];
598 for(j=0; (j<ny); j++) {
599 tmp[YY] = j*box[YY];
600 for(k=0; (k<nz); k++) {
601 tmp[ZZ] = k*box[ZZ];
602 for(m=0; (m<natom); m++,n++) {
603 if ((!bOrdered && ((m % natmol) == 0)) || bOrdered)
604 rvec_add(xx[n % natom],tmp,xx[n]);
605 else
606 ;
607 }
608 }
609 }
610 }
611
612 clear_mat(boxje);
613 boxje[XX][XX] = box[XX]*nx;
614 boxje[YY][YY] = box[YY]*ny;
615 boxje[ZZ][ZZ] = box[ZZ]*nz;
616
617 printf("Crystal: %10.5f %10.5f %10.5f\n",
618 nx*box[XX],ny*box[YY],nz*box[ZZ]);
619
620 if (debug && !bDiamond) {
621 if (bSeries)
622 for(i=3; (i<=10*rcut); i++) {
623 fprintf(debug,"This is with rcut = %g\n",i*0.1);
624 virial(debug,bFull,nmax/natmol,xx,boxje,
625 0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
626 }
627 else
628 virial(debug,bFull,nmax/natmol,xx,boxje,
629 rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
630 }
631
632 if (bDiamond)
633 mk_diamond(pdba,xx,odist,&symtab,bPBC,boxje);
634
635 fn = ftp2fn(efSTO,NFILE,fnm);
636 if (fn2ftp(fn) == efPDB) {
637 fp = gmx_ffopen(fn,"w");
638 if (bDiamond)
639 fprintf(fp,"HEADER This is a *diamond*\n");
640 else
641 fprintf(fp,"HEADER A beautiful Ice Ih crystal\n");
642 fprintf(fp,"REMARK Generated by mkice with the following options:\n"
643 "REMARK nx = %d, ny = %d, nz = %d, odist = %g, hdist = %g\n",
644 nx,ny,nz,odist,hdist);
645 bromacs(quote,255);
646 write_pdbfile(fp,quote,pdba,xx,boxje,' ',-1);
647 gmx_ffclose(fp);
648 }
649 else {
650 bromacs(quote,255);
651 write_sto_conf(fn,quote,pdba,xx,NULL,boxje);
652 }
653
654 if (ftp2bSet(efTRN,NFILE,fnm))
655 write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
656
657 return 0;
658 }
659