|
24 | 24 | # |
25 | 25 | #-- |
26 | 26 |
|
27 | | -from argparse import ArgumentParser |
28 | | -import os |
29 | | - |
30 | | -import h5py as h5 |
31 | | - |
32 | | -from molmod.io.chk import load_chk |
33 | | -from yaff import System |
34 | | - |
35 | | -from quickff.tools import set_ffatypes, get_ei_radii, average, charges_to_bcis |
36 | | -from quickff.io import read_abinitio, make_yaff_ei, read_bci_constraints |
37 | | - |
38 | | - |
39 | | -description = '''\ |
40 | | -This script reads atomic charges from an input file and makes a Yaff parameters file |
41 | | -suitable for the QuickFF option --ei.''' |
42 | | - |
43 | | - |
44 | | -def parse_args(args=None): |
45 | | - parser = ArgumentParser(description=description) |
46 | | - parser.add_argument( |
47 | | - '-v', '--verbose', default=False, action='store_true', |
48 | | - help='Increase verbosity of the script [default=%(default)s].' |
49 | | - ) |
50 | | - parser.add_argument( |
51 | | - '--ffatypes', default=None, |
52 | | - choices=['None','list_of_atypes', 'low','medium','high','highest'], |
53 | | - help='Assign atom types in the system by parsing an ordered list of' |
54 | | - 'atom types as argument or through the automatic built-in ' |
55 | | - 'detection (see documentation). By default (or if None is given), ' |
56 | | - 'the atom types are assumed to be defined in the input files. ' |
57 | | - '[default=%(default)s]' |
58 | | - ) |
59 | | - parser.add_argument( |
60 | | - '--gaussian', default=False, action='store_true', |
61 | | - help='Use gaussian smeared charges. The radii are taken from the input ' |
62 | | - 'file fn_in (from dataset /path/radii for HDF5 file or from label ' |
63 | | - '`radii` for CHK file) if the data is present, otherwise the radii ' |
64 | | - 'are estimated according to the procedure of Chen et al. See ' |
65 | | - '``quickff.tools.get_ei_radii`` for more info.' |
66 | | - ) |
67 | | - parser.add_argument( |
68 | | - '--bci', default=False, action='store_true', |
69 | | - help='Convert averaged atomic charges to bond charge increments, i.e. ' |
70 | | - 'charge transfers along the chemical bonds in the system. In this ' |
71 | | - 'way, one is certain of ending up with a globally neutral system ' |
72 | | - 'even if bcis from different systems are combined. This option ' |
73 | | - 'requires the definition of bonds, if the bonds cannot be read ' |
74 | | - 'from the system file, they will be estimated from the interatomic ' |
75 | | - 'distances using the detect_bonds routine in the Yaff System ' |
76 | | - 'class. [default=%(default)s]' |
77 | | - ) |
78 | | - parser.add_argument( |
79 | | - '--bci-constraints', default=None, |
80 | | - help='A file containing constraints for the charge to bci fit in a ' |
81 | | - 'master: slave0,slave1,...: sign format. A new line should be used ' |
82 | | - 'for each master and the format is insensitive towards spaces.' |
83 | | - 'Sign should be 1.0 or -1.0 indicating wheter or not a sign switch ' |
84 | | - 'should be introduced when mapping the slaves to the master.' |
85 | | - ) |
86 | | - parser.add_argument( |
87 | | - 'fn_sys', |
88 | | - help='Any file from which the system can be extracted (MolMod CHK, Gaussian ' |
89 | | - 'FCHK, XYZ, ...).' |
90 | | - ) |
91 | | - parser.add_argument( |
92 | | - 'charges', |
93 | | - help='The atomic charges to be used. This argument has the form fn_charges:path, ' |
94 | | - 'where fn_charges is a filename that contains the charges and path refers ' |
95 | | - 'to a location within that file where the charges can be found. If ' |
96 | | - 'fn_charges is an HDF5 file, path is the location of a dataset containing ' |
97 | | - 'the charges, e.g. \'/charges\'. If fn_charges is a MolMod CHK file, path ' |
98 | | - 'is the label of the dataset that contains the atomic charges.' |
99 | | - ) |
100 | | - parser.add_argument( |
101 | | - 'fn_out', default='pars_ei.txt', nargs='?', |
102 | | - help='Name of the Yaff file to write the parameters to. [default=%(default)s]' |
103 | | - ) |
104 | | - if args is None: |
105 | | - args = parser.parse_args() |
106 | | - else: |
107 | | - args = parser.parse_args(args.split()) |
108 | | - if args.charges.count(':') != 1: |
109 | | - parser.error('The argument charges must contain exactly one colon.') |
110 | | - return args |
111 | | - |
112 | | - |
113 | | -def main(args=None): |
114 | | - if args is None: |
115 | | - args = parse_args() |
116 | | - else: |
117 | | - args = parse_args(args) |
118 | | - # Load system file |
119 | | - if args.fn_sys.endswith('.fchk'): |
120 | | - numbers, coords, energy, grad, hess, masses, rvecs, pbc = read_abinitio(args.fn_sys, do_hess=False) |
121 | | - system = System(numbers, coords, rvecs=None, charges=None, radii=None, masses=masses) |
122 | | - system.detect_bonds() |
123 | | - else: |
124 | | - system = System.from_file(args.fn_sys) |
125 | | - |
126 | | - # Guess atom types if needed |
127 | | - if args.ffatypes is not None: |
128 | | - set_ffatypes(system, args.ffatypes) |
129 | | - ffatypes = [system.ffatypes[i] for i in system.ffatype_ids] |
130 | | - |
131 | | - # Load atomic charges |
132 | | - fn_charges, _, path = args.charges.partition(':') |
133 | | - if fn_charges.endswith('.h5'): |
134 | | - with h5.File(fn_charges, 'r') as f: |
135 | | - if not path in f: |
136 | | - raise IOError('Given HDF5 file %s does not contain a dataset %s' % (fn_charges, path)) |
137 | | - charges = f[path][:] |
138 | | - radii = None |
139 | | - if args.gaussian: |
140 | | - path_radii = os.path.join(os.path.dirname(path), 'radii') |
141 | | - if 'radii' in f[path]: |
142 | | - radii = average(f['%s/radii' %path][:], ffatypes, fmt='dict') |
143 | | - else: |
144 | | - radii = average(get_ei_radii(system.numbers), ffatypes, fmt='dict') |
145 | | - elif fn_charges.endswith('.chk'): |
146 | | - sample = load_chk(fn_charges) |
147 | | - if path in list(sample.keys()): |
148 | | - charges = sample[path] |
149 | | - else: |
150 | | - raise IOError('Given CHK file %s does not contain a dataset with label %s' % (fn_charges, path)) |
151 | | - radii = None |
152 | | - if args.gaussian: |
153 | | - if 'radii' in list(sample.keys()): |
154 | | - radii = average(sample['radii'], ffatypes, fmt='dict') |
155 | | - else: |
156 | | - raise IOError('Invalid extension, fn_charges should be a HDF5 or a CHK file.') |
157 | | - |
158 | | - # Derive charge parameters |
159 | | - if args.bci: |
160 | | - constraints = {} |
161 | | - if args.bci_constraints is not None: |
162 | | - constraints = read_bci_constraints(args.bci_constraints) |
163 | | - bcis = charges_to_bcis(charges, ffatypes, system.bonds, constraints=constraints, verbose=args.verbose) |
164 | | - make_yaff_ei(args.fn_out, None, bcis=bcis, radii=radii) |
165 | | - else: |
166 | | - charges = average(charges, ffatypes, fmt='dict', verbose=args.verbose) |
167 | | - make_yaff_ei(args.fn_out, charges, radii=radii) |
168 | | - |
| 27 | +from quickff.scripts import qff_input_ei |
169 | 28 |
|
170 | 29 | if __name__=='__main__': |
171 | | - main() |
| 30 | + qff_input_ei() |
0 commit comments