Demo entry 6660255

asdf

   

Submitted by anonymous on Nov 14, 2017 at 16:57
Language: Python. Code size: 12.0 kB.

#!/usr/bin/env desres-exec
#{
# exec desres-cleanenv -c --user-env        \
#   -m Python/2.7.1-06A/bin                 \
#   -m destro/0.6.2/lib-python              \
#   -m scipy/0.12.0-01/lib-python           \
#   -m destiny/0.4.47-3/lib-python          \
#   -m molpro/2012.1.12-desres-3-02/serial/bin \
#   -m psi4/4.0b5-desres-2-01/bin           \
#   -- python -u $0 "$@"
#}

import sys
sys.path.append('/u/en/donchev/pyprogs')
sys.path.append('/u/en/donchev/pyprogs/chemsys')
sys.path.append('/d/en/ffdevel-0/bin/pyqmdata')
sys.path.append('/u/en/donchev/playground/qmdata6/')

import math, numpy
import subprocess
import time
import copy
import string

import units
import prepmpr
import pth4 as pth
import frame
import class_qmjob
import class_scannmer



class MultiScanNmer(object):
    def __init__(self, finput, opts, molpro_args=[]):
        self.molpro_args = molpro_args
        self.opts = opts
        #self.basopt = self.opts.basopt
        #self.basene = self.opts.basene
        #self.basesl = self.opts.basesl
        self.ncores = self.opts.ncores
        self.verbose = self.opts.verbose
        self.rerun  = self.opts.rerun
        self.pinput = pth.Path(finput)
        self.folder = self.pinput.folder
        self.folder.cd()
        self.name = self.pinput.name()
        self.setup = self.pinput.eval()
        self.progress_file = self.pinput.chext('.results')
        self.progress_lines = []
        self.dict_file = self.pinput.chext('.dict')
        #self.myz_file = self.pinput.chext('.myz')
        #self.selection_file = self.pinput.chext('.selection')

        self.start_time = time.time()

        self.pydict = pth.Path(self.setup['dict']).eval()
        self.frames = {}
        for k,v in self.pydict['data'].items():#[:1]:
            self.frames[k] = frame.LightMolecule(elements=self.pydict['elements'], xyzxyz=v['xyz'], energy=v['relenergy'])
            self.frames[k].split_monomers([int(x) for x in self.pydict['natoms'].split()], [int(x) for x in self.pydict['charges'].split()])
            self.frames[k].k = k
            #moncnfs = [int(x) for x in v['moncnfs'].split('-')]
            self.frames[k].data = v
            self.frames[k].scan = class_scannmer.ScanNmer(self.frames[k], self.folder+'scan'+str(k), self.opts)

        #_v = [(v.energy,k,v) for k,v in self.frames.items()]
        #self.frames = {min(_v)[1]:min(_v)[2], max(_v)[1]:max(_v)[2]}
        #print min(_v)[:2], max(_v)[:2]

        self.error_file = self.folder + 'error.py'
        self.error_archive_file = self.folder + 'errarchive.py'


    ################################################################################################
    def run(self):
        jobs = sum([f.scan.jobs for f in self.frames.values()],[])
        info = self.jobs2info(jobs, rerun=self.rerun > 0)

        for f in self.frames.values():
            f.scan.run()

        if self.opts.correct_solv:
            for f in self.frames.values():
                f.scan.run_field()
            #return
            jobs = sum([f.scan.get_jobs_solvated_correction() for f in self.frames.values()],[])
            info = self.jobs2info(jobs, rerun=self.rerun > 0)
            #print info
            for f in self.frames.values():
                f.scan.run_solvated_correction()

        self.save_conformers()
        self.rewrite_progress()


    ################################################################################################
    def save_conformers(self):
        info = copy.copy(self.pydict)
        minene = min(f.scan.minene for f in self.frames.values())
        minmonene = min(f.scan.monenergy for f in self.frames.values())

        info['data'] = {}
        data_lines = []
        for k,f in sorted(self.frames.items()):
            s = f.scan
            scan_info_fmt = "'relenergy':{:.5f}, 'intenergy':{:.5f}, 'trans_amides':{:5s}, 'lnpop':{:8.5f}, 'monenergy':{:.5f}, 'moncnfs':'{}'"
            scan_info = scan_info_fmt.format(s.minene-minene, s.minintene, str(f.data['trans_amides']), f.data['lnpop'], s.monenergy-minmonene, f.data['moncnfs'])
            if self.opts.solvated:
                scan_info += ", 'solute':{}".format(s.SOLUTE)
            data_lines.append("%3d: {" % k + scan_info)
            for idsf,sf in sorted(s.scan.items()):
                data_lines.append(self.one_line_dict_info(idsf,sf))
            data_lines[-1] += '}'

        txt = self.folder.pprint_dict_ordered(info, head='title  nmols  natoms  charges  effsize'.split(), tail='species  maeff'.split())
        txt = txt.replace('{}', '{\n'+',\n'.join(data_lines)+'}')
        self.dict_file.write(txt)


    #################################################################################################
    def one_line_xyz_info(self, info, f):
        fmt = '%.3f inte=%.3f nat=%s fch=%s tam=%s eff=%d mids=%s'
        vals = (f.energy, f.data['intenergy'], ':'.join([x for x in info['natoms'].split()]), ':'.join([x for x in info['charges'].split()]), f.data['trans_amides'], info['effsize'], f.data['moncnfs'])
        return fmt % vals


    #################################################################################################
    def one_line_dict_info(self, i, f):
        def p(s): return f.ec[s] * units.au2kcal #/ self.pydict['nmols']
        if self.opts.solvated:
            fmt = "{:3d}: dict(solv={:12.5e}, solv2body={:12.5e}, solvnonadd={:12.5e}, dimer_int={}, xyz='{}')"
        else:
            fmt = "{:3d}: dict(inter={:12.5e}, e2body={:12.5e}, nonadd={:12.5e}, dimer_int={}, xyz='{}')"
        return fmt.format(i, p('eintr'), p('ene2body'), p('enonadd'), str(f.ec['prettyintedimer']), f.xyzxyz)


    ################################################################################################
    def jobs2info(self, jobs, rerun=True, data_mode=None):
        class_qmjob.process_in_parallel(copy.copy(jobs), 'run()', self.ncores, delay=0.)

        errors = [j.error for j in jobs if j.error is not None]
        if errors:
            self.save_errors(errors)
            print 'found errors', str(self.pinput)
            assert not errors, errors

        self.error_file.rm()
        return dict([(j.label,j.data) for j in jobs])


    ################################################################################################
    @property
    def error_archive(self):
        if not self.error_archive_file.exists: return []
        try:
            return self.error_archive_file.eval()
        except SyntaxError:
            return []


    ################################################################################################
    def save_errors(self, errors):
        self.error_file.rm()
        self.error_file.pprint([], errors)
        errarch = self.error_archive + errors
        self.error_archive_file.pprint([], errarch)


    ################################################################################################
    def rewrite_progress(self):
        lines  = ['molpro basis sets: small = %s  large = %s' % (self.opts.basissmall, self.opts.basislarge)]
        lines += ['ghost options: rcut_ghost = %.2f  fast_nonadd = %s' % (self.opts.rcut_ghost, self.opts.fast_nonadd)]
        lines += ['']
        for k,f in sorted(self.frames.items()):
            lines += f.scan.progress_lines
        self.progress_file.write(lines)


    ################################################################################################
    def clean_after_yourself(self):
        fld = self.pinput.folder
        for f in fld.glob('*.mpr') + fld.glob('*.out*') + fld.glob('*.xml*') + fld.glob('*log*'):
            f.rm()





if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Monomer minimizer with (H)XH constraints")
    parser.add_argument('input_py', action='store',
        help='Input file with a Python dictionary setting')
    parser.add_argument('-n', '--ncores', action='store', default=4, type=int,
        help='Number of cores for parallel execution')
    parser.add_argument('--preclean', action='store_true', default=False,
        help='Delete all intermidiate molpro files BEFORE running QM')
    parser.add_argument('--clean', action='store_true', default=False,
        help='Delete all intermidiate molpro files AFTER running QM')
    parser.add_argument('--basissmall', action='store', choices='vdz avdz vtz avtz vqz avqz'.split(), default='avdz',
        help='Single-point energy basis set for optimized geometry, default %(default)s; will be converted to desX-rev5')
    parser.add_argument('--basislarge', action='store', choices='vtz avtz vqz avqz v5z av5z'.split(), default='avqz',
        help='Single-point energy basis set for optimized geometry, default %(default)s; will be converted to desX-rev5')
    parser.add_argument('--rerun', action='store', default=0, type=int,
        help='If 0, try to avoid reruning QM')
    parser.add_argument('-r', '--rcut_ghost', action='store', default=100., type=float,
        help='Cutoff radius for placing ghost atom when computing dimers in the small basis set; if the radius is small and --fast_nonadd option is set to True, interaction and nonadditive energies might be contaminated with BSSE (multimer beyond dimers)')
    parser.add_argument('--fast_nonadd', action='store_true', default=False,
        help='Try placing as few ghost atoms as possible when computing dimer energies in the small basis set; interaction and nonadditive energies contaminated with BSSE (multimer beyond dimers)')
    parser.add_argument('--schedule', action='store', default='range(-5,6)', type=str,
        help='Schedule for the breathing scan; in units of ~0.1 Anstrom')
    parser.add_argument('-v', '--verbose', action='store_true', default=False,
        help='Report progress on screen')
    parser.add_argument('--solvated', action='store_true', default=False,
        help='Compute solvation energy rather than interaction energy; assume that all waters comprise the solvent')
    parser.add_argument('--correct_solv', action='store_true', default=False,
        help='Use bigger basis to compute nonadditive component of solvation energy for the cluster of waters experiencing strong electric field exerted by the solute; assume that all waters comprise the solvent')
    parser.add_argument('--basismedium', action='store', choices='vdz avdz vtz avtz vqz avqz'.split(), default='avtz',
        help='Single-point energy basis set for optimized geometry, default %(default)s; will be converted to desX-rev5')
    parser.add_argument('--minwatc', action='store', default=6, type=int,
        help='Use at least this many waters to estimate the medium basis correction to nonadditive energy')
    parser.add_argument('--maxwatc', action='store', default=10, type=int,
        help='Use no more than this many waters to estimate the medium basis correction to nonadditive energy')
    parser.add_argument('--field_thresh', action='store', default=0.02, type=float,
        help='Use all waters experiencing at least this strong electric field to estimate the medium basis correction to nonadditive energy')

    args = parser.parse_args()
    args.schedule = eval(args.schedule)
    #print args
    if args.correct_solv:
        assert args.basismedium != args.basissmall

    #molpro_args = []
    #if args.d is not None: molpro_args += ['-d',args.d]
    #if args.I is not None: molpro_args += ['-I',args.I]
    #if args.W is not None: molpro_args += ['-W',args.W]
    #if args.L is not None: molpro_args += ['-L',args.L]

    job = MultiScanNmer(args.input_py, args)
    if args.preclean: job.clean_after_yourself()
    job.run()
    if args.clean: job.clean_after_yourself()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).