OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
waterSphere
1#!@Python3_EXECUTABLE@
2
3__author__ = "Chris Fennell and Dan Gezelter"
4__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
5__license__ = "OpenMD"
6
7import sys
8import os
9import math
10import random
11import argparse
12from argparse import RawDescriptionHelpFormatter
13
14def usage():
15 print(__doc__)
16
17# function to calculate the arccosine
18def acos(rad):
19 return math.atan2(math.sqrt(1 - rad * rad), rad)
20
21# Function to check if a point is within the spherical boundaries
22def inSphere(x, y, z, radius):
23 return math.sqrt(x**2 + y**2 + z**2) <= radius
24
25# Function to open the output file
26def open_output_file(fileName):
27 f = open(fileName, 'w')
28 if os.path.exists(fileName):
29 return f
30 else:
31 sys.exit("Error: can't open file {}\n".format(fileName))
32
33# Function to find the cutoff
34def find_cutoff(cutoff, boxx, boxy, boxz):
35 bm = min(boxx, boxy, boxz)
36 boxLength2 = 0.5 * bm
37 if boxLength2 <= cutoff:
38 cutoff = int(boxLength2)
39
40# Functions to print the water model types
41def print_cl(outfile):
42 outfile.write("""
43molecule{
44 name = "Cl-";
45 atom[0]{
46 type = "Cl-";
47 position(0.0, 0.0, 0.0);
48 }
49}""")
50
51def print_na(outfile):
52 outfile.write("""
53molecule{
54 name = "Na+";
55 atom[0]{
56 type = "Na+";
57 position(0.0, 0.0, 0.0);
58 }
59}""")
60
61def print_ssd_e(outfile):
62 outfile.write("""
63molecule{
64 name = "SSD_E";
65 atom[0]{
66 type = "SSD_E";
67 position( 0.0, 0.0, 0.0 );
68 orientation( 0.0, 0.0, 0.0 );
69 }
70}""")
71
72def print_ssd_rf(outfile):
73 outfile.write("""
74molecule{
75 name = "SSD_RF";
76 atom[0]{
77 type = "SSD_RF";
78 position( 0.0, 0.0, 0.0 );
79 orientation( 0.0, 0.0, 0.0 );
80 }
81}""")
82
83def print_ssd(outfile):
84 outfile.write("""
85molecule{
86 name = "SSD";
87 atom[0]{
88 type = "SSD";
89 position( 0.0, 0.0, 0.0 );
90 orientation( 0.0, 0.0, 0.0 );
91 }
92}""")
93
94def print_ssdq(outfile):
95 outfile.write("""
96molecule{
97 name = "SSDQ";
98 atom[0]{
99 type = "SSDQ";
100 position( 0.0, 0.0, 0.0 );
101 orientation( 0.0, 0.0, 0.0 );
102 }
103}""")
104
105def print_ssdqo(outfile):
106 outfile.write("""
107molecule{
108 name = "SSDQO";
109 atom[0]{
110 type = "SSDQO";
111 position( 0.0, 0.0, 0.0 );
112 orientation( 0.0, 0.0, 0.0 );
113 }
114}""")
115
116def print_ssd1(outfile):
117 outfile.write("""
118molecule{
119 name = "SSD1";
120 atom[0]{
121 type = "SSD1";
122 position( 0.0, 0.0, 0.0 );
123 orientation( 0.0, 0.0, 0.0 );
124 }
125}""")
126
127def print_tip3p(outfile):
128 outfile.write("""
129molecule{
130 name = "TIP3P";
131 atom[0]{
132 type = "O_TIP3P";
133 position( 0.0, 0.0, -0.06556 );
134 }
135 atom[1]{
136 type = "H_TIP3P";
137 position( 0.0, 0.75695, 0.52032 );
138 }
139 atom[2]{
140 type = "H_TIP3P";
141 position( 0.0, -0.75695, 0.52032 );
142 }
143 rigidBody[0]{
144 members(0, 1, 2);
145 }
146}""")
147
148def print_tip4p(outfile):
149 outfile.write("""
150molecule{
151 name = "TIP4P";
152 atom[0]{
153 type = "O_TIP4P";
154 position( 0.0, 0.0, -0.06556 );
155 }
156 atom[1]{
157 type = "H_TIP4P";
158 position( 0.0, 0.75695, 0.52032 );
159 }
160 atom[2]{
161 type = "H_TIP4P";
162 position( 0.0, -0.75695, 0.52032 );
163 }
164 atom[3]{
165 type = "EP_TIP4P";
166 position( 0.0, 0.0, 0.08444 );
167 }
168 rigidBody[0]{
169 members(0, 1, 2, 3);
170 }
171}""")
172
173def print_tip4p_ice(outfile):
174 outfile.write("""
175molecule{
176 name = "TIP4P-Ice";
177 atom[0]{
178 type = "O_TIP4P-Ice";
179 position( 0.0, 0.0, -0.06556 );
180 }
181 atom[1]{
182 type = "H_TIP4P-Ice";
183 position( 0.0, 0.75695, 0.52032 );
184 }
185 atom[2]{
186 type = "H_TIP4P-Ice";
187 position( 0.0, -0.75695, 0.52032 );
188 }
189 atom[3]{
190 type = "EP_TIP4P-Ice";
191 position( 0.0, 0.0, 0.09214 );
192 }
193 rigidBody[0]{
194 members(0, 1, 2, 3);
195 }
196}""")
197
198def print_tip4p_2005(outfile):
199 outfile.write("""
200molecule{
201 name = "TIP4P-2005";
202 atom[0]{
203 type = "O_TIP4P-2005";
204 position( 0.0, 0.0, -0.06556 );
205 }
206 atom[1]{
207 type = "H_TIP4P-2005";
208 position( 0.0, 0.75695, 0.52032 );
209 }
210 atom[2]{
211 type = "H_TIP4P-2005";
212 position( 0.0, -0.75695, 0.52032 );
213 }
214 atom[3]{
215 type = "EP_TIP4P-2005";
216 position( 0.0, 0.0, 0.08904 );
217 }
218 rigidBody[0]{
219 members(0, 1, 2, 3);
220 }
221}""")
222
223def print_tip4pew(outfile):
224 outfile.write("""
225molecule{
226 name = "TIP4P-Ew";
227 atom[0]{
228 type = "O_TIP4P-Ew";
229 position( 0.0, 0.0, -0.06556 );
230 }
231 atom[1]{
232 type = "H_TIP4P-Ew";
233 position( 0.0, 0.75695, 0.52032 );
234 }
235 atom[2]{
236 type = "H_TIP4P-Ew";
237 position( 0.0, -0.75695, 0.52032 );
238 }
239 atom[3]{
240 type = "EP_TIP4P-Ew";
241 position( 0.0, 0.0, 0.05944 );
242 }
243 rigidBody[0]{
244 members(0, 1, 2, 3);
245 }
246}""")
247
248def print_tip5p(outfile):
249 outfile.write("""
250molecule{
251 name = "TIP5P";
252 atom[0]{
253 type = "O_TIP5P";
254 position( 0.0, 0.0, -0.06556 );
255 }
256 atom[1]{
257 type = "H_TIP5P";
258 position( 0.0, 0.75695, 0.52032 );
259 }
260 atom[2]{
261 type = "H_TIP5P";
262 position( 0.0, -0.75695, 0.52032 );
263 }
264 atom[3]{
265 type = "EP_TIP5P";
266 position( 0.57154, 0.0, -0.46971 );
267 }
268 atom[4]{
269 type = "EP_TIP5P";
270 position( -0.57154, 0.0, -0.46971 );
271 }
272 rigidBody[0]{
273 members(0, 1, 2, 3, 4);
274 }
275}""")
276
277def print_tip5pe(outfile):
278 outfile.write("""
279molecule{
280 name = "TIP5P-E";
281 atom[0]{
282 type = "O_TIP5P-E";
283 position( 0.0, 0.0, -0.06556 );
284 }
285 atom[1]{
286 type = "H_TIP5P";
287 position( 0.0, 0.75695, 0.52032 );
288 }
289 atom[2]{
290 type = "H_TIP5P";
291 position( 0.0, -0.75695, 0.52032 );
292 }
293 atom[3]{
294 type = "EP_TIP5P";
295 position( 0.57154, 0.0, -0.46971 );
296 }
297 atom[4]{
298 type = "EP_TIP5P";
299 position( -0.57154, 0.0, -0.46971 );
300 }
301 rigidBody[0]{
302 members(0, 1, 2, 3, 4);
303 }
304}""")
305
306def print_spce(outfile):
307 outfile.write("""
308molecule{
309 name = "SPCE";
310 atom[0]{
311 type = "O_SPCE";
312 position( 0.0, 0.0, -0.06461 );
313 }
314 atom[1]{
315 type = "H_SPCE";
316 position( 0.0, 0.81649, 0.51275 );
317 }
318 atom[2]{
319 type = "H_SPCE";
320 position( 0.0, -0.81649, 0.51275 );
321 }
322 rigidBody[0]{
323 members(0, 1, 2);
324 }
325}""")
326
327def print_spc(outfile):
328 outfile.write("""
329molecule{
330 name = "SPC";
331 atom[0]{
332 type = "O_SPC";
333 position( 0.0, 0.0, -0.06461 );
334 }
335 atom[1]{
336 type = "H_SPC";
337 position( 0.0, 0.81649, 0.51275 );
338 }
339 atom[2]{
340 type = "H_SPC";
341 position( 0.0, -0.81649, 0.51275 );
342 }
343 rigidBody[0]{
344 members(0, 1, 2);
345 }
346}""")
347
348def print_spc_hw(outfile):
349 outfile.write("""
350molecule{
351 name = "SPC-HW";
352 atom[0]{
353 type = "O_SPC-HW";
354 position( 0.0, 0.0, -0.06461 );
355 }
356 atom[1]{
357 type = "D_SPC-HW";
358 position( 0.0, 0.81649, 0.51275 );
359 }
360 atom[2]{
361 type = "D_SPC-HW";
362 position( 0.0, -0.81649, 0.51275 );
363 }
364 rigidBody[0]{
365 members(0, 1, 2);
366 }
367}""")
368
369def print_ne6(outfile):
370 outfile.write("""
371molecule{
372 name = "NE6";
373 atom[0]{
374 type = "O_NE6";
375 position( 0.0, 0.0, 0.0 );
376 }
377 atom[1]{
378 type = "H_NE6";
379 position( 0.0, 0.576029, 0.79283665 );
380 }
381 atom[2]{
382 type = "H_NE6";
383 position( 0.0, -0.576029, 0.79283665 );
384 }
385 atom[3]{
386 type = "EP_NE6";
387 position( 0.0, 0.23, 0.0 );
388 }
389 atom[4]{
390 type = "LP_NE6";
391 position( 0.732813007, -0.50364843, 0.0 );
392 }
393 atom[5]{
394 type = "LP_NE6";
395 position( -0.732813007, -0.50364843, 0.0 );
396 }
397 rigidBody[0]{
398 members(0, 1, 2, 3, 4, 5);
399 }
400}""")
401
402def print_tip3p_fb(outfile):
403 outfile.write("""
404molecule{
405 name = "TIP3P-FB";
406 atom[0]{
407 type = "O_TIP3P-FB";
408 position( 0.0, 0.0, -0.066424 );
409 }
410 atom[1]{
411 type = "H_TIP3P-FB";
412 position( 0.0, 0.819341, 0.527225 );
413 }
414 atom[2]{
415 type = "H_TIP3P-FB";
416 position( 0.0, -0.819341, 0.527225 );
417 }
418 rigidBody[0]{
419 members(0, 1, 2);
420 }
421}""")
422
423def print_tip4p_fb(outfile):
424 outfile.write("""
425molecule{
426 name = "TIP4P-FB";
427 atom[0]{
428 type = "O_TIP4P-FB";
429 position( 0.0, 0.0, -0.0655549 );
430 }
431 atom[1]{
432 type = "H_TIP4P-FB";
433 position( 0.0, 0.75695, 0.520327 );
434 }
435 atom[2]{
436 type = "H_TIP4P-FB";
437 position( 0.0, -0.75695, 0.520327 );
438 }
439 atom[3]{
440 type = "EP_TIP4P-FB";
441 position( 0.0, 0.0, 0.0397151 );
442 }
443 rigidBody[0]{
444 members(0, 1, 2, 3);
445 }
446}""")
447
448def print_opc(outfile):
449 outfile.write("""
450molecule{
451 name = "OPC";
452 atom[0]{
453 type = "O_OPC";
454 position( 0.0, 0.0, -0.0603651 );
455 }
456 atom[1]{
457 type = "H_OPC";
458 position( 0.0, 0.685582, 0.479134 );
459 }
460 atom[2]{
461 type = "H_OPC";
462 position( 0.0, -0.685582, 0.479134 );
463 }
464 atom[3]{
465 type = "EP_OPC";
466 position( 0.0, 0.0, 0.0990349 );
467 }
468 rigidBody[0]{
469 members(0, 1, 2, 3);
470 }
471}""")
472
473def print_opc3(outfile):
474 outfile.write("""
475molecule{
476 name = "OPC3";
477 atom[0]{
478 type = "O_OPC3";
479 position( 0.0, 0.0, -0.0632382 );
480 }
481 atom[1]{
482 type = "H_OPC3";
483 position( 0.0, 0.799262, 0.501939 );
484 }
485 atom[2]{
486 type = "H_OPC3";
487 position( 0.0, -0.799262, 0.501939 );
488 }
489 rigidBody[0]{
490 members(0, 1, 2);
491 }
492}""")
493
494def print_dpd(outfile):
495 outfile.write("""
496molecule{
497 name = "DPD";
498 atom[0]{
499 type = "DPD";
500 position( 0.0, 0.0, 0.0 );
501 }
502}""")
503
504def print_cg2(outfile):
505 outfile.write("""
506molecule{
507 name = "CG2";
508 atom[0]{
509 type = "CG2";
510 position( 0.0, 0.0, 0.0 );
511 }
512}""")
513
514# Function to print fake water model
515def print_fake_water(waterName, outfile):
516 outfile.write("\n\nmolecule{{\n name = \"{}\";\n atom[0]{{\n type = \"{}\";\n position(0.0, 0.0, 0.0);\n }}\n}}".format(waterName, waterName))
517
518# Function to print the water model
519def print_water_model(waterCase, outfile):
520 if waterCase == 0: print_cl(outfile)
521 elif waterCase == 1: print_na(outfile)
522 elif waterCase == 2: print_ssd_e(outfile)
523 elif waterCase == 3: print_ssd_rf(outfile)
524 elif waterCase == 4: print_ssd(outfile)
525 elif waterCase == 5: print_ssd1(outfile)
526 elif waterCase == 6: print_tip3p(outfile)
527 elif waterCase == 7: print_tip4p(outfile)
528 elif waterCase == 8: print_tip4pew(outfile)
529 elif waterCase == 9: print_tip5p(outfile)
530 elif waterCase == 10: print_tip5pe(outfile)
531 elif waterCase == 11: print_spce(outfile)
532 elif waterCase == 12: print_spc(outfile)
533 elif waterCase == 13: print_dpd(outfile)
534 elif waterCase == 14: print_cg2(outfile)
535 elif waterCase == 15: print_ssdq(outfile)
536 elif waterCase == 16: print_ssdqo(outfile)
537 elif waterCase == 17: print_tip4p_ice(outfile)
538 elif waterCase == 18: print_tip4p_2005(outfile)
539 elif waterCase == 19: print_spc_hw(outfile)
540 elif waterCase == 20: print_ne6(outfile)
541 elif waterCase == 21: print_tip3p_fb(outfile)
542 elif waterCase == 22: print_tip4p_fb(outfile)
543 elif waterCase == 23: print_opc(outfile)
544 elif waterCase == 24: print_opc3(outfile)
545
546# Function to validate the water model
547def validate_water(waterName):
548 water_models = {
549 'Cl-': 0, 'Na+': 1, 'SSD_E': 2, 'SSD_RF': 3, 'SSD': 4, 'SSD1': 5,
550 'TIP3P': 6, 'TIP4P': 7, 'TIP4P-Ew': 8, 'TIP5P': 9, 'TIP5P-E': 10,
551 'SPCE': 11, 'SPC': 12, 'DPD': 13, 'CG2': 14, 'SSDQ': 15, 'SSDQO': 16,
552 'TIP4P-Ice': 17, 'TIP4P-2005': 18, 'SPC-HW': 19, 'NE6': 20,
553 'TIP3P-FB': 21, 'TIP4P-FB': 22, 'OPC': 23, 'OPC3': 24
554 }
555 waterCase = water_models.get(waterName, -1)
556 invalidWater = False
557 if waterCase == -1:
558 invalidWater = True
559 return (waterCase, invalidWater)
560
561# Function to print the meta data
562def print_meta_data(outfile, doWaterInc, waterName, waterCase, invalidWater,
563 nMol, cutoff, alpha):
564 outfile.write(" <MetaData>\n")
565
566 if doWaterInc:
567 outfile.write('#include "water.inc"')
568 else:
569 print_water_model(waterCase, outfile)
570
571 if invalidWater:
572 print_fake_water(waterName, outfile)
573
574 outfile.write("\n\n")
575 outfile.write("component{\n")
576 outfile.write(" type = \"{}\";\n".format(waterName))
577 outfile.write(" nMol = {};\n".format(nMol))
578 outfile.write("}\n")
579 outfile.write("ensemble = \"LangevinHull\";\n")
580 outfile.write("forceField = \"Water\";\n")
581 outfile.write("electrostaticSummationMethod = \"shifted_force\";\n")
582 outfile.write("electrostaticScreeningMethod = \"damped\";\n")
583 outfile.write("cutoffRadius = {};\n".format(cutoff))
584 outfile.write("dampingAlpha = {};\n".format(alpha))
585 outfile.write("usePeriodicBoundaryConditions = \"false\";\n")
586 outfile.write("\n")
587 outfile.write("targetTemp = 300;\n")
588 outfile.write("targetPressure = 1.0;\n")
589 outfile.write("\n")
590 outfile.write("tauThermostat = 1e3;\n")
591 outfile.write("tauBarostat = 1e4;\n")
592 outfile.write("\n")
593 outfile.write("dt = 2.0;\n")
594 outfile.write("runTime = 1e3;\n")
595 outfile.write("\n")
596 outfile.write("tempSet = \"true\";\n")
597 outfile.write("thermalTime = 10;\n")
598 outfile.write("sampleTime = 100;\n")
599 outfile.write("statusTime = 2;\n")
600 outfile.write(" </MetaData>\n")
601
602# Function to print the frame data
603def print_frame_data(outfile, boxx, boxy, boxz):
604 outfile.write(" <Snapshot>\n")
605 outfile.write(" <FrameData>\n")
606 outfile.write(" Time: 0\n")
607 outfile.write(" Hmat: {{{{ {}, 0, 0 }}, {{ 0, {}, 0 }}, {{ 0, 0, {} }}}}\n".format(boxx, boxy, boxz))
608 outfile.write(" </FrameData>\n")
609
610# Function to write the output file
611def write_output_file(outfile, nMol, xCorr, yCorr, zCorr, cutoff, alpha,
612 boxx, boxy, boxz, waterName, waterCase, invalidWater,
613 doRandomize, doWaterInc):
614 # Write out the header
615 outfile.write("<OpenMD version=2>\n")
616 find_cutoff(cutoff, boxx, boxy, boxz)
617 print_meta_data(outfile, doWaterInc, waterName, waterCase, invalidWater,
618 nMol, cutoff, alpha)
619 print_frame_data(outfile, boxx, boxy, boxz)
620 outfile.write(" <StuntDoubles>\n")
621
622 # Shift the box center to the origin and write out the coordinates
623 for i in range(nMol):
624 xCorr[i] -= 0.5 * boxx
625 yCorr[i] -= 0.5 * boxy
626 zCorr[i] -= 0.5 * boxz
627
628 q0 = 1.0
629 q1 = q2 = q3 = 0.0
630
631 if doRandomize:
632 cosTheta = 2.0 * random.random() - 1.0
633 theta = acos(cosTheta)
634 phi = 2.0 * math.pi * random.random()
635 psi = 2.0 * math.pi * random.random()
636
637 q0 = math.cos(0.5 * theta) * math.cos(0.5 * (phi + psi))
638 q1 = math.sin(0.5 * theta) * math.cos(0.5 * (phi - psi))
639 q2 = math.sin(0.5 * theta) * math.sin(0.5 * (phi - psi))
640 q3 = math.cos(0.5 * theta) * math.sin(0.5 * (phi + psi))
641
642 outfile.write("{}\tpq\t{} {} {} {} {} {} {}\n".format(i, xCorr[i], yCorr[i], zCorr[i], q0, q1, q2, q3))
643
644 outfile.write(" </StuntDoubles>\n")
645 outfile.write(" </Snapshot>\n")
646 outfile.write("</OpenMD>\n")
647
648def main(argv):
649
650 parser = argparse.ArgumentParser(
651 description='builds spheres of water',
652 formatter_class=argparse.RawDescriptionHelpFormatter,
653 epilog=
654'''
655Note: you can only use values of -O that are smaller
656 than the derived radius for a given density and
657 number of molecules.
658
659Example:
660 waterSphere -d 0.997 -I 20 -O 40 -w SSD_RF -o ssdrfWater.omd
661''')
662 parser.add_argument("-d", "--density=", action="store", dest='density',
663 type=float, help="density in g/cm^3",
664 required=False, default=1.0)
665 parser.add_argument("-l", "--lattice=", choices=[0,1],
666 action="store", dest='lattice',
667 type=int, default=0, required=False,
668 help="0 - face centered cubic, 1 - simple cubic")
669 parser.add_argument("-c", "--rcut=", action="store", dest="rcut",
670 type=float, default=3.2, required=False,
671 help="default overlap cutoff in angstroms")
672 parser.add_argument("-o", "--output=", action="store", dest='fileName',
673 help="output file name", required=False,
674 default="freshWater.omd")
675 parser.add_argument("-w", "--water=", action="store", dest="waterName",
676 help="name of the water StuntDouble",
677 required=False, default="SPCE")
678 parser.add_argument("-r", "--randomize",
679 action=argparse.BooleanOptionalAction, default=False,
680 type=bool, dest='doRandomize',
681 help="randomize orientations")
682 parser.add_argument("-O", "--outer_radius", action="store", required=True,
683 dest='outer_radius', type=float,
684 help="outer radius of the water sphere")
685 parser.add_argument("-I", "--inner_radius", action="store", required=False,
686 dest='inner_radius', type=float, default=0.0,
687 help="inner radius of the water sphere")
688 parser.add_argument("-m", "--printInclude",
689 action=argparse.BooleanOptionalAction, default=False,
690 type=bool, dest='doWaterInc',
691 help="print out a water.inc file (file with all water models)")
692
693 if len(sys.argv) == 1:
694 parser.print_help()
695 sys.exit(2)
696 args = parser.parse_args()
697
698 density = args.density
699 lattice = args.lattice
700 rcut = args.rcut
701 fileName = args.fileName
702 waterName = args.waterName
703 outer_radius = args.outer_radius
704 inner_radius = args.inner_radius
705 boxx = boxy = boxz = 2 * outer_radius
706 doRandomize = args.doRandomize
707 doWaterInc = args.doWaterInc
708
709 tolerance = 1.0E-8
710 mass = 2.99151E-23 # mass of H2O in grams
711 cm3ToAng3 = 1E24 # convert cm^3 to angstroms^3
712 densityConvert = mass * cm3ToAng3
713 cutoff = 9
714 alpha = 0.18
715 invalidWater = False
716 waterCase = -1
717
718 if doWaterInc:
719 if fileName == 'water.inc':
720 sys.exit("Error: {} cannot be \"water.inc\"\n Please choose a different name\n".format(fileName))
721 if density <= 0:
722 usage()
723 sys.exit("Error: the value for '-d' ({}) is not a valid number\n Please choose a positive real # value\n".format(density))
724 if rcut <= 0:
725 usage()
726 sys.exit("Error: the value for '-c' ({}) is not a valid number\n Please choose a positive real # value\n".format(arg))
727 if lattice != 0 and lattice != 1:
728 usage()
729 sys.exit("Error: the '-l' value ({}) is not a valid number\n Please choose 0 or 1\n".format(arg))
730 if outer_radius <= 0:
731 usage()
732 sys.exit("Error: the value for '-O' ({}) is not a valid number\n Please choose a positive real # value\n".format(outer_radius))
733 if inner_radius < 0:
734 usage()
735 sys.exit("Error: the value for '-I' ({}) is not a valid number\n Please choose a positive real # value\n".format(inner_radius))
736
737 outfile = open_output_file(fileName)
738
739 # Set magic lattice numbers
740 if lattice == 0:
741 a = (4 * 18.01 / (0.602 * density)) ** (1.0 / 3.0)
742 acut = rcut * math.sqrt(2.0)
743 if acut > a:
744 a = acut
745 newDensity = 4.0 * 18.01 / (0.602 * a ** 3)
746 print("using density of {} to match cutoff value ({})\n".format(newDensity, rcut))
747 elif lattice == 1:
748 a = (18.01 / (0.602 * density)) ** (1.0 / 3.0)
749 acut = rcut
750 if acut > a:
751 a = acut
752 newDensity = 18.01 / (0.602 * a ** 3)
753 print("using density of {} to match cutoff value ({})\n".format(newDensity, rcut))
754
755 nxFloat = boxx / a
756 nx = int(nxFloat + nxFloat / abs(nxFloat * 2))
757 nyFloat = boxy / a
758 ny = int(nyFloat + nyFloat / abs(nyFloat * 2))
759 nzFloat = boxz / a
760 nz = int(nzFloat + nzFloat / abs(nzFloat * 2))
761
762 anew = min(boxx / nx, boxy / ny, boxz / nz)
763
764 if anew < acut:
765 anew = acut
766
767 nxFloat = boxx / anew
768 nx = int(nxFloat + nxFloat / abs(nxFloat * 2))
769 nyFloat = boxy / anew
770 ny = int(nyFloat + nyFloat / abs(nyFloat * 2))
771 nzFloat = boxz / anew
772 nz = int(nzFloat + nzFloat / abs(nzFloat * 2))
773
774 if lattice == 0:
775 nMol = 4 * nx * ny * nz
776 else:
777 nMol = nx * ny * nz
778
779 newDensity = nMol * densityConvert / (boxx * boxy * boxz)
780
781 if abs(newDensity - density) > tolerance:
782 print("Resetting density to {} to make chosen box sides work out\n".format(newDensity))
783
784 cellLengthX = boxx / nx
785 cellLengthY = boxy / ny
786 cellLengthZ = boxz / nz
787
788 cell2X = cellLengthX * 0.5
789 cell2Y = cellLengthY * 0.5
790 cell2Z = cellLengthZ * 0.5
791
792 xCorr = []
793 yCorr = []
794 zCorr = []
795
796 centerx = 0.5 * boxx
797 centery = 0.5 * boxy
798 centerz = 0.5 * boxz
799
800
801 if lattice == 0:
802 xref = [0.0, 0.0, cell2X, cell2X]
803 yref = [0.0, cell2Y, cell2Y, 0.0]
804 zref = [0.0, cell2Z, 0.0, cell2Z]
805
806 counter = 0
807 for z in range(nz):
808 for y in range(ny):
809 for x in range(nx):
810 for uc in range(4):
811 xt = xref[uc] + cellLengthX * x - centerx
812 yt = yref[uc] + cellLengthY * y - centery
813 zt = zref[uc] + cellLengthZ * z - centerz
814 if inSphere(xt, yt, zt, outer_radius) and \
815 inSphere(xt, yt, zt, inner_radius) == False:
816 xCorr.append(xt)
817 yCorr.append(yt)
818 zCorr.append(zt)
819 counter += 1
820 elif lattice == 1:
821 xref = [cell2X]
822 yref = [cell2Y]
823 zref = [cell2Z]
824
825 counter = 0
826 for z in range(nz):
827 for y in range(ny):
828 for x in range(nx):
829 xt = xref[uc] + cellLengthX * x - centerx
830 yt = yref[uc] + cellLengthY * y - centery
831 zt = zref[uc] + cellLengthZ * z - centerz
832 if inSphere(xt, yt, zt, outer_radius) and \
833 inSphere(xt, yt, zt, inner_radius) == False:
834 xCorr.append(x_val)
835 yCorr.append(y_val)
836 zCorr.append(z_val)
837 counter += 1
838
839 nMol = counter
840
841 (waterCase, invalidWater) = validate_water(waterName)
842 if invalidWater:
843 print("Warning: '{}' is not a recognized water model name.\n".format(waterName))
844 print(" Use the '-m' option to generate a 'water.inc' with the\n")
845 print(" recognized water model geometries.\n\n")
846 if waterName == 'DPD':
847 density *= 0.25
848 elif waterName == 'CG2':
849 density *= 0.5
850
851 write_output_file(outfile, nMol, xCorr, yCorr, zCorr, cutoff, alpha,
852 boxx, boxy, boxz, waterName, waterCase, invalidWater,
853 doRandomize, doWaterInc)
854 outfile.close()
855
856 print("The water box \"{}\" was generated.\n".format(fileName))
857
858 if doWaterInc:
859 try:
860 outfile = open("./water.inc", "w")
861 except:
862 sys.exit("Error: can't open file {}\n".format("water.inc"))
863
864 outfile.write("#ifndef _WATER_INC_\n#define _WATER_INC_\n")
865 print_cl(outfile)
866 print_na(outfile)
867 print_ssd_e(outfile)
868 print_ssd_rf(outfile)
869 print_ssd(outfile)
870 print_ssd1(outfile)
871 print_tip3p(outfile)
872 print_tip4p(outfile)
873 print_tip4pew(outfile)
874 print_tip5p(outfile)
875 print_tip5pe(outfile)
876 print_spce(outfile)
877 print_spc(outfile)
878 print_dpd(outfile)
879 print_cg2(outfile)
880 print_ssdq(outfile)
881 print_ssdqo(outfile)
882 print_tip4p_ice(outfile)
883 print_tip4p_2005(outfile)
884 print_spc_hw(outfile)
885 print_ne6(outfile)
886 print_tip3p_fb(outfile)
887 print_tip4p_fb(outfile)
888 print_opc(outfile)
889 print_opc3(outfile)
890 outfile.write("\n\n#endif")
891
892 print("The file \"water.inc\" was generated for inclusion in \"{}\"\n".format(fileName))
893
894 outfile.close()
895
896if __name__ == "__main__":
897 main(sys.argv[1:])