MRPyCM/computeSolvent.py at main · Gabrielgerez/MRPyCM · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
from vampyr import vampyr3d as vp
import MRPyCM as mpcm
import numpy as np


if __name__ == '__main__':
    import sys


# Global parameters
def run(*args, **kwargs):
    """
    input parameter dictionary:
    "order" : int
    "box" : list of floats
    "prec" : float
    "charge_width" : float
    "charges" : list of floats
    "charge_coords" : list of lists of floats
    "cav_coords" : list of lists of floats
    "cav_radii" : list of floats
    "boundary_width" : float
    "eps_out" : float
    "perm_formulation" : string
    "solver_type" : string
    "ionic_strength" : float
    "max_iter" : int
    "kain_hist" : int
    """
    # Define parameters and defaults
    keys =kwargs.keys()
    k = kwargs["order"] if ("order" in keys) else 5                            # Polynomial order
    L = kwargs["box"] if ("box" in keys) else  [-10,10]                        # Simulation box size
    epsilon = kwargs["prec"] if ("prec" in keys) else 1.0e-4                   # Relative precision

    charge_width = kwargs["charge_width"] if ("charge_width" in keys) else 1000.0
    charges = kwargs["charges"] if ("charges" in keys) else [1.0]
    charge_coords = kwargs["charge_coords"] if ("charge_coords" in keys) else [[0.0000000000,    0.0000000000,    0.000000000]]

    cav_coords = kwargs["cav_coords"] if ("cav_coords" in keys) else  charge_coords
    cav_radii = kwargs["cav_radii"] if ("cav_radii" in keys) else [1.0]
    boundary_width = kwargs["boundary_width"] if ("boundary_width" in keys) else 0.2

    eps_out = kwargs["eps_out"] if ("eps_out" in keys) else 2.0
    perm_formulation = kwargs["perm_formulation"] if "perm_formulation" in keys else "exponential"

    solvent_type = kwargs["solver_type"] if ("solver_type" in keys) else "gpe"
    ionic_strength = kwargs["ionic_strength"] if ("ionic_strength" in keys) else 0.1

    max_iter = kwargs["max_iter"] if ("max_iter" in keys) else 100
    kain_hist = kwargs["kain_hist"] if ("kain_hist" in keys) else 0


    # Define MRA and multiwavelet projector
    MRA = vp.MultiResolutionAnalysis(order=k, box=L)
    print(MRA)
    P_eps = vp.ScalingProjector(mra=MRA, prec=epsilon)
    D_abgv = vp.ABGVDerivative(mra=MRA, a=0.0, b=0.0)
    Poissop = vp.PoissonOperator(mra=MRA, prec=epsilon)

    # nuclear density and total molecular density to compute the vacuum potential
    dens = P_eps(mpcm.constructChargeDensity(charge_coords, charges, width_parameter=charge_width))

    # Solvent part
    C = mpcm.Cavity(cav_coords, cav_radii, boundary_width)

    if ("linear" == perm_formulation.lower()):
        perm = P_eps(mpcm.LinPerm(C, inside=1.0, outside=eps_out))
    else:
        perm = P_eps(mpcm.ExpPerm(C, inside=1.0, outside=eps_out))

    if ("pb" == solvent_type.lower()):
        k_sq = P_eps(mpcm.DHScreening(C, inside=0.0, outside=mpcm.computeKappaOut(eps_out, ionic_strength)))
        Solver = mpcm.PBSolver(dens, perm, k_sq,  Poissop, D_abgv, max_iter=max_iter, hist=kain_hist)

    elif ("lpb" == solvent_type.lower()):
        k_sq = P_eps(mpcm.DHScreening(C, inside=0.0, outside=mpcm.computeKappaOut(eps_out, ionic_strength)))
        Solver = mpcm.LPBSolver(dens, perm, k_sq, Poissop, D_abgv, max_iter=max_iter, hist=kain_hist)

    else:
        Solver = mpcm.GPESolver(dens, perm, Poissop, D_abgv, max_iter=max_iter, hist=kain_hist)


    reaction_op = mpcm.ReactionOperator(Solver)
    reaction_op.setup(epsilon)
    E_R = reaction_op.trace()
    print("E_R: ", E_R)
    return E_R, Solver.iterations


if __name__ == '__main__':
    arg_dict = eval(sys.argv[1])
    energy, iterations = run(**arg_dict)
    print("Energy: ", energy, " Iterations: ", iterations)