@@ -1,15 +1,10 @@
Loading
1 -
2 1
from respyte.molecule import *
3 2
from respyte.readinp_resp import Input
4 3
from respyte.molecule_resp import *
5 4
from respyte.select_grid import *
6 -
7 5
from respyte.resp_unit import *
8 6
9 -
10 7
def main():
11 -
12 -
    # cwd = current working directory in which input folder exists
13 8
    cwd = os.getcwd()
14 9
    # check if the current working idr contains input folder
15 10
    print('\n')
@@ -17,7 +12,7 @@
Loading
17 12
    print(' 1. Reading  input files and folders.                    ')
18 13
    print('---------------------------------------------------------')
19 14
    if os.path.isdir("%s/input" % cwd):
20 -
        print(' Found the input folder. Now read input/input.yml')
15 +
        print(' Found the input folder. Now read input/respyte.yml')
21 16
    else:
22 17
        print(' Failed to find input folder. Should have input(folder) containing input.yml and molecules(folder)')
23 18
    # read respyte.yml
@@ -31,6 +26,7 @@
Loading
31 26
    else:
32 27
        print(' Not using Cheminformatics?')
33 28
        molecule = Molecule_respyte()
29 +
34 30
    molecule.addInp(inp)
35 31
36 32
    for idx, i in enumerate(inp.nmols):
@@ -77,20 +73,16 @@
Loading
77 73
                        if (len(numbers)==4):
78 74
                            xyz = [x for x in numbers[0:3]]
79 75
                            pts.append(xyz)
80 -
                selectedPts = SelectGridPts(tmpmol,inner,outer,pts, radii )
81 -
                print(' Read %d pts from %s_%s.espf' % (len(pts), molN, confN))
82 -
                print(' select pts: inner =%0.4f, outer = %0.4f' %(inner, outer))
83 -
                print(' Use %d pts out of %d pts' %(len(selectedPts), len(pts)))
76 +
                selectedPtsIdx, selectedPts = SelectGridPts(tmpmol,inner,outer,pts, radii )
84 77
            else:
85 -
                selectedPts = None
86 -
            listofselectedPts.append(selectedPts)
78 +
                selectedPtsIdx = None
79 +
            listofselectedPts.append(selectedPtsIdx)
87 80
        molecule.addCoordFiles(*coordfilepath)
88 81
        molecule.addEspf(*espffilepath, selectedPts = listofselectedPts)
89 82
    print('---------------------------------------------------------')
90 83
    print(' 2. Charge fitting to QM data                            ')
91 84
    print('---------------------------------------------------------')
92 85
    os.chdir(cwd)
93 -
    print(' resp calculation is running on %s' % cwd)
94 86
    cal = Respyte_Optimizer()
95 87
    cal.addInp(inp)
96 88
    cal.addMolecule(molecule)
@@ -102,19 +94,29 @@
Loading
102 94
            aval = inp.restraintinfo['a']
103 95
            bval = inp.restraintinfo['b']
104 96
            cal.Model3qpot(aval, bval)
97 +
        elif inp.restraintinfo['penalty'] == 'model4':
98 +
            a1val = inp.restraintinfo['a1']
99 +
            a2val = inp.restraintinfo['a2']
100 +
            bval = inp.restraintinfo['b']
101 +
            cal.Model4qpot(a1val, a2val, bval)
102 +
        elif inp.restraintinfo['penalty'] == 'model5':
103 +
            a1val = inp.restraintinfo['a1']
104 +
            a2val = inp.restraintinfo['a2']
105 +
            bval = inp.restraintinfo['b']
106 +
            cal.Model5qpot(a1val, a2val, bval)
107 +
        elif inp.restraintinfo['penalty'] == 'model6':
108 +
            aval = inp.restraintinfo['a']
109 +
            bval = inp.restraintinfo['b']
110 +
            cal.Model6qpot(aval, bval)
105 111
        elif inp.restraintinfo['penalty'] == '2-stg-fit':
106 112
            a1val = inp.restraintinfo['a1']
107 113
            a2val = inp.restraintinfo['a2']
108 114
            bval = inp.restraintinfo['b']
109 -
            cal.twoStageFit(a1val, a2val,bval)
110 -
111 -
    # print()
112 -
    # print('    #####################################################')
113 -
    # print('    ###               Test calculations               ###')
114 -
    # print('    #####################################################')
115 -
    # cal.Model2qpot(0.005)
116 -
    # cal.Model3qpot(0.0005,0.1)
117 -
    # cal.twoStageFit(0.0005,0.001,0.1)
115 +
            cal.twoStageFit(a1val, a2val, bval)
116 +
        elif inp.restraintinfo['penalty'] == '2-stg-fit(6)':
117 +
            aval = inp.restraintinfo['a']
118 +
            bval = inp.restraintinfo['b']
119 +
            cal.twoStageFit_Model6(aval, bval)
118 120
119 121
if __name__ == '__main__':
120 122
    main()

@@ -143,10 +143,10 @@
Loading
143 143
            # copy espf file generated in moli/confj/tmp to moli/confj/
144 144
            copyfile(espfoutput, path + '/%s_%s.espf' % (molN, confN))
145 145
            copyfile(griddat, path + '/grid.dat')
146 -
            print('\n')
147 -
            print('---------------------------------------------------------')
148 -
            print('       Done! Hope it helps, Will see you again:)         ')
149 -
            print('---------------------------------------------------------')
146 +
    print('\n')
147 +
    print('---------------------------------------------------------')
148 +
    print('       Done! Hope it helps, Will see you again:)         ')
149 +
    print('---------------------------------------------------------')
150 150
151 151
if __name__ == "__main__":
152 152
    main()

@@ -37,6 +37,7 @@
Loading
37 37
        self.resnumbers = []
38 38
        #self.equivGroups = [] # maybe don't need this
39 39
        self.listofpolars = []
40 +
        self.listofburieds = []
40 41
        self.listofchargeinfo = []
41 42
42 43
        self.gridxyzs = []
@@ -89,7 +90,10 @@
Loading
89 90
    def addInp(self, inpcls):
90 91
        assert isinstance(inpcls, Input)
91 92
        self.inp = inpcls
92 -
93 +
    def changesymmetry(self, symmetry=False):
94 +
        print(self.inp.symmetry)
95 +
        self.inp.symmetry = symmetry
96 +
        print('After:',self.inp.symmetry)
93 97
    def removeSingleElemList(self, lst):
94 98
        needtoremove = []
95 99
        for idx, i in enumerate(lst):
@@ -298,6 +302,8 @@
Loading
298 302
            # For now, when cheminformatics is not used, ignore polar atoms
299 303
            listofpolar = []
300 304
            self.listofpolars.append(listofpolar)
305 +
            listofburied = []
306 +
            self.listofburieds.append(listofburied)
301 307
302 308
    def addEspf(self, *espfFiles, selectedPts):
303 309
        for idx, espfFile in enumerate(espfFiles):
@@ -362,6 +368,7 @@
Loading
362 368
            self.atomnames.append([])
363 369
            self.resnumbers.append([])
364 370
            self.listofpolars.append([])
371 +
            self.listofburieds.append([])
365 372
            xyzs = []
366 373
            self.nmols.append(len(xyzs))
367 374
            indices = []
@@ -438,11 +445,15 @@
Loading
438 445
                    # Using openeye tool, make listofpolar,
439 446
                    symmetryClass = []
440 447
                    listofpolar = []
448 +
                    listofburied = []
441 449
                    oechem.OEAssignHybridization(oemol)
442 450
                    for atom in oemol.GetAtoms():
443 451
                        symmetryClass.append(atom.GetSymmetryClass())
444 452
                        if atom.IsCarbon() and int(atom.GetHyb()) != 3:
453 +
445 454
                            listofpolar.append(atom.GetIdx())
455 +
                        if len([bond for bond in atom.GetBonds()]) >3:
456 +
                            listofburied.append(atom.GetIdx())
446 457
                            # ispolar = False
447 458
                            # for bond in atom.GetBonds():
448 459
                            #     atom2 = bond.GetNbr(atom)
@@ -461,16 +472,13 @@
Loading
461 472
                                    listofpolar.append(atom.GetIdx())
462 473
                                elif atom2.IsCarbon() and atom2.GetIdx() in listofpolar:
463 474
                                    listofpolar.append(atom.GetIdx())
475 +
                                if atom2.GetIdx() in listofburied:
476 +
                                    listofburied.append(atom.GetIdx()) # store hydrogens bonded to buried atoms
464 477
                        if atom.IsPolar():
465 478
                            listofpolar.append(atom.GetIdx())
466 479
467 -
                            # for atom2 in oemol.GetAtoms():
468 -
                            #     if atom2.IsHydrogen() and atom2.IsConnected(atom) is True:
469 -
                            #         listofpolar.append(atom2.GetIdx())
470 -
                            #     elif atom2.IsCarbon() and atom2.IsConnected(atom) and oechem.OEGetHybridization(atom2) < 3:
471 -
                            #         listofpolar.append(atom2.GetIdx())
472 -
473 480
                    listofpolar = sorted(set(listofpolar))
481 +
                    listofburied = sorted(set(listofburied))
474 482
                    idxof1statm, resnameof1statm = self.getidxof1statm(resnumber, resname)
475 483
                    unique_resid = set(resnameof1statm)
476 484
                    sameresid = [[i for i, v in enumerate(resnameof1statm) if v == value] for value in unique_resid]
@@ -527,13 +535,27 @@
Loading
527 535
                            for j in self.atomidinfo[i]:
528 536
                                newatomidinfo[newid].append(j)
529 537
                            del newatomidinfo[i]
530 -
                    self.atomids.append(newatomid)
531 -
                    self.atomidinfo = newatomidinfo
538 +
                    # print('newatomidinfo', newatomidinfo)
539 +
                    # print('oldatomidinfo', self.atomidinfo)
540 +
                    # print('newatomid', newatomid)
541 +
                    # print('oldatomid', atomid)
542 +
                    if self.inp is not None:
543 +
                        print(self.inp.symmetry)
544 +
                        if self.inp.symmetry == False:
545 +
                            self.atomids.append(atomid)
546 +
                            self.atomidinfo = self.atomidinfo
547 +
                        else:
548 +
                            self.atomids.append(newatomid)
549 +
                            self.atomidinfo = newatomidinfo
550 +
                    else:
551 +
                        self.atomids.append(newatomid)
552 +
                        self.atomidinfo = newatomidinfo
532 553
                    self.elems.append(atomicNum)
533 554
                    self.resnames.append(resname)
534 555
                    self.atomnames.append(atomname)
535 556
                    self.resnumbers.append(resnumber)
536 557
                    self.listofpolars.append(listofpolar)
558 +
                    self.listofburieds.append(listofburied)
537 559
538 560
                    if self.inp is not None:
539 561
                        chargeinfo = self.gen_chargeinfo(self.inp.resChargeDict, newatomid, self.atomidinfo, resnumber)
@@ -710,8 +732,16 @@
Loading
710 732
                            for j in self.atomidinfo[i]:
711 733
                                newatomidinfo[newid].append(j)
712 734
                            del newatomidinfo[i]
713 -
                    self.atomids.append(newatomid)
714 -
                    self.atomidinfo = newatomidinfo
735 +
                    if self.inp is not None:
736 +
                        if self.inp.symmetry == False:
737 +
                            self.atomids.append(atomid)
738 +
                            self.atomidinfo = self.atomidinfo
739 +
                        else:
740 +
                            self.atomids.append(newatomid)
741 +
                            self.atomidinfo = newatomidinfo
742 +
                    else:
743 +
                        self.atomids.append(newatomid)
744 +
                        self.atomidinfo = newatomidinfo
715 745
                    self.elems.append(atomicNum)
716 746
                    self.resnames.append(resname)
717 747
                    self.atomnames.append(atomname)
@@ -788,8 +818,8 @@
Loading
788 818
        # print('atom ids:    ', molecule.atomids,'\n')
789 819
        # #print('resnames:    ',molecule.resnames,'\n')
790 820
        # print('polar atoms: ', molecule.listofpolars[-1],'\n')
791 -
        #print('charge info: ',molecule.listofchargeinfo[-1],'\n')
792 -
        #print('atomidinfo', molecule.atomidinfo)
821 +
        # #print('charge info: ',molecule.listofchargeinfo[-1],'\n')
822 +
        # #print('atomidinfo', molecule.atomidinfo)
793 823
        #print('-----------------------------------------------------------------------------------------')
794 824
    # print(molecule.nmols, len(molecule.nmols))
795 825
    # print(len(molecule.elems), len(molecule.atomids), len(molecule.resnames), len(molecule.listofpolars), len(molecule.xyzs))

@@ -1,5 +1,4 @@
Loading
1 1
import os,sys, copy
2 -
import scipy.stats as stats
3 2
import scipy as sci
4 3
import numpy as np
5 4
import re
@@ -13,7 +12,6 @@
Loading
13 12
    import openeye.oechem as oechem
14 13
except ImportError:
15 14
    warn(' The Openeye module cannot be imported. ( Please provide equivGoups and listofpolar manually.)')
16 -
17 15
from respyte.molecule import *
18 16
from respyte.readinp_resp import Input
19 17
@@ -21,17 +19,18 @@
Loading
21 19
bohr2Ang = 0.52918825 # change unit from bohr to angstrom
22 20
23 21
class Respyte_Optimizer:
24 -
    def __init__(self, inp = None, molecule = None):
22 +
    def __init__(self, inp = None, molecule = None, normalization = False):
25 23
        self.inp = inp
26 24
        self.molecule = molecule
27 -
25 +
        self.normalization = normalization
28 26
    def addInp(self, inpcls):
29 27
        assert isinstance(inpcls, Input)
30 28
        self.inp = inpcls
31 -
29 +
        self.normalization = inpcls.normalization
32 30
    def addMolecule(self, molcls):
33 31
        self.molecule = molcls
34 -
32 +
    def Normalization(self, normalization):
33 +
        self.normalization = normalization
35 34
    def EspDesignMatrix(self, xyzs, gridxyzs, espval):
36 35
        """
37 36
        Produces a design matrix A and vector B using esp values to solve for charges q in Aq=B.
@@ -99,6 +98,79 @@
Loading
99 98
100 99
        return apot, bpot
101 100
101 +
    def NewEspDesignMatrix(self, xyzs, gridxyzs, espval, efval, objN):
102 +
        nAtoms = xyzs.shape[0]
103 +
        molbohr = np.array(xyzs.copy())
104 +
        gridx = np.array([xyz[0] for xyz in gridxyzs])
105 +
        gridy = np.array([xyz[1] for xyz in gridxyzs])
106 +
        gridz = np.array([xyz[2] for xyz in gridxyzs])
107 +
        potv = np.array(espval)
108 +
        potvsq = potv*potv
109 +
        ssvpot = np.sum( potvsq)
110 +
        ssvpotsorted = np.sum( np.sort(potvsq) )
111 +
112 +
        abs_ef = []
113 +
        squared_ef = []
114 +
        for vec in efval:
115 +
            absef = np.sqrt(np.dot(np.array(vec), np.array(vec)))
116 +
            abs_ef.append(absef)
117 +
            squared_ef.append(absef*absef)
118 +
        abs_ef = np.array(abs_ef)
119 +
        squared_ef = np.array(squared_ef)
120 +
        squared_esp = potvsq
121 +
        abs_esp = np.sqrt(squared_esp)
122 +
123 +
        if objN == 'obj2':
124 +
            weights = abs_esp
125 +
        elif objN == 'obj2_2':
126 +
            weights = (abs_esp**2 + 0.01**2 )**(1/2)
127 +
        elif objN == 'obj3':
128 +
            weights = squared_esp
129 +
        elif objN == 'obj3_2':
130 +
            weights = (squared_esp**2 + 0.01**2 )**(1/2)
131 +
        elif objN == 'obj4':
132 +
            weights = abs_ef
133 +
        elif objN == 'obj4_2':
134 +
            weights = (abs_ef**2 + 0.05**2 )**(1/2)
135 +
        elif objN == 'obj5':
136 +
            weights = squared_ef
137 +
        else:
138 +
            raise NotImplementedError('%s is not implemented. Try obj2, obj3, obj4, or obj5 instead!' % objN)
139 +
        invRijSq = []
140 +
        ainvRijSq = []
141 +
        invRij = []
142 +
        ainvRij = []
143 +
        dxij = []
144 +
        dyij = []
145 +
        dzij = []
146 +
        for atom in range(nAtoms):
147 +
            idx = atom
148 +
            dxi = gridx-molbohr[idx,0]
149 +
            dxij.append( dxi)
150 +
            dyi = gridy-molbohr[idx,1]
151 +
            dyij.append( dyi)
152 +
            dzi = gridz-molbohr[idx,2]
153 +
            dzij.append( dzi)
154 +
            rijsq = dxi*dxi + dyi*dyi +dzi*dzi
155 +
            invRiSq = 1.0/rijsq
156 +
            invRijSq.append( invRiSq )
157 +
            ainvRijSq.append(weights* invRiSq)
158 +
            invRij.append( np.sqrt( invRiSq) )
159 +
            ainvRij.append(weights*np.sqrt( invRiSq))
160 +
        # build A matrix and B vector
161 +
        apot = np.zeros( (nAtoms, nAtoms) )
162 +
        bpot = np.zeros( nAtoms)
163 +
        for j, ainvRi in enumerate( ainvRij):
164 +
            bi = potv*ainvRi
165 +
            bpot[j] = np.sum( np.sort(bi) )
166 +
            for k in range(0, j):
167 +
                sumrijrik = np.sum( np.sort(ainvRi*invRij[k]) )
168 +
                apot[j,k] = sumrijrik
169 +
                apot[k,j] = sumrijrik
170 +
        for j, ainvRiSq in enumerate( ainvRijSq):
171 +
            apot[j,j] = np.sum(ainvRiSq)
172 +
173 +
        return apot, bpot
102 174
    def EfDesignMatrix(self, xyzs, gridxyzs, efval):
103 175
        """
104 176
        Produces a design matrix A and vector B from electric field values to solve for charges q in Aq=B.
@@ -181,17 +253,9 @@
Loading
181 253
        fldv = np.array(efval)
182 254
        fldvsq = fldv*fldv
183 255
        ssvfld = np.sum( fldvsq)
184 -
        # weight = np.sqrt(ssvpot/ssvfld)
185 -
        # print('weight:',ssvpot/ssvfld,', sqrt of weight:',np.sqrt(ssvpot/ssvfld))
186 -
        # print('compare As:',len(Aesp), Aesp[3],len(Aef), Aef[3])
187 -
        # print('compare Bs:',len(Besp), Besp, len(Bef), Bef)
188 256
        weight = ssvpot/ssvfld
189 -
        Acomb =( Aesp + Aef*weight)
190 -
        Bcomb =( Besp + Bef*weight)
191 -
        # print('A after mixin:', Acomb[3])
192 -
        # print('B after mixin:', Bcomb)
193 -
        # input()
194 -
257 +
        Acomb =( Aesp + Aef*weight)/2
258 +
        Bcomb =( Besp + Bef*weight)/2
195 259
        return Acomb, Bcomb
196 260
197 261
    def LagrangeChargeConstraint(self, aInp, bInp, chargeinfo):
@@ -221,7 +285,6 @@
Loading
221 285
        charges = []
222 286
        check = []
223 287
        unique_resname = list(set(chargeinfo[i][1] for i,info in enumerate(chargeinfo)))
224 -
        # print(unique_resname)
225 288
        for lstofidx, resname, netchg in chargeinfo:
226 289
            if len(lstofidx) == 0 :
227 290
                unique_resname.remove(resname)
@@ -245,8 +308,6 @@
Loading
245 308
                for i in range(nAtom, nAtom+len(unique_resname)):  # Need to remove duplicates
246 309
                    lastrow.append(0.0)
247 310
                lastrows.append(lastrow)
248 -
        # print(list(len(lastrows[i]) for i in range(len(lastrows))))
249 -
        # print('Here', aInp.shape)
250 311
        apot = []
251 312
        for i, row in enumerate(aInp):
252 313
            newrow = list(row)
@@ -352,10 +413,7 @@
Loading
352 413
        crossProd = np.sum(np.dot(qpot, bpot))
353 414
        modelSumSq = np.dot(qpot,np.dot(apot,qpot))
354 415
        chiSq = sumSq - 2*crossProd + modelSumSq
355 -
        print(' max(Aij) : ', np.amax(apot))
356 -
        print(' chi esp  : ', chiSq)
357 -
        scalingFac = (float(self.inp.gridspace)/0.7)**3
358 -
        print(' chi esp after Norm : ', chiSq*scalingFac)
416 +
        chiSq_rstr = np.sum((qpot**2 + (0.1)**2)**(0.5)-0.1)
359 417
        rrmsval = np.sqrt(chiSq/sumSq)
360 418
        return rrmsval
361 419
@@ -384,7 +442,6 @@
Loading
384 442
        apot = copy.deepcopy(np.array(apot))
385 443
        bpot = copy.deepcopy(np.array(bpot))
386 444
        qpot = copy.deepcopy(np.array(qpot))
387 -
        # sumSq = np.sum(np.dot(espval,espval))
388 445
        fldv = np.array(efval)
389 446
        fldvsq = fldv*fldv
390 447
        ssvfld = np.sum( fldvsq)
@@ -394,9 +451,43 @@
Loading
394 451
        rrmsval = np.sqrt(chiSq/ssvfld)
395 452
        return rrmsval
396 453
397 -
##########################################################################################################
398 -
###               Functions for avoiding any duplicates in different fitting models                    ###
399 -
##########################################################################################################
454 +
    def centerofmass(self, xyzs, elems):
455 +
        M = sum(float(list(PeriodicTable.values())[int(elems[i]-1)]) for i in range(len(elems)))
456 +
        return np.sum([xyzs[i,:] *float(list(PeriodicTable.values())[int(elems[i]-1)]) / M for i in range(len(xyzs))],axis=0)
457 +
458 +
    def MMdipole(self, xyzs, elems, qpot):
459 +
        com = self.centerofmass(xyzs, elems)
460 +
        mux = 0
461 +
        muy = 0
462 +
        muz = 0
463 +
        for qi, xyzi in zip(qpot, xyzs):
464 +
            mux += qi*(xyzi[0]-com[0])
465 +
            muy += qi*(xyzi[1]-com[1])
466 +
            muz += qi*(xyzi[2]-com[2])
467 +
            ri_squared = np.dot(np.array(xyzi)-com,np.array(xyzi)-com)
468 +
        mu_squared = mux**2+ muy**2+muz**2
469 +
        mu =  np.sqrt(mu_squared)*bohr2ang/0.20819434
470 +
        return mu
471 +
472 +
    def Ncond(self, apotInp, elemInp, qpotInp, Options):
473 +
        AA = copy.deepcopy(apotInp)
474 +
        for i in range(len(elemInp)):
475 +
            if elemInp[i] == 'H' or elemInp[i] == 1:
476 +
                continue
477 +
            else:
478 +
                if 'weights' in Options and 'tightness' in Options:
479 +
                    weights = Options['weights']
480 +
                    AA[i][i] -= weights[i]*qpotInp[i] / np.sqrt(qpotInp[i]**2 + float(Options['tightness'])**2)**3
481 +
                else:
482 +
                    continue
483 +
        try:
484 +
            w, v = np.linalg.eig(AA)
485 +
            cond = max(abs(w))/min(abs(w))
486 +
        except:
487 +
            print('condition number calculation failed.')
488 +
            cond = 0
489 +
        return cond
490 +
400 491
    def combineMatrices(self):
401 492
        loc = 0
402 493
        apots = []
@@ -404,37 +495,50 @@
Loading
404 495
        atomid_comb = []
405 496
        elem_comb = []
406 497
        Size = 0
498 +
        add = 0
407 499
        for idx, i in enumerate(self.molecule.nmols):
408 500
            if i == 0:
409 -
                pass
501 +
                add += 1
410 502
            else:
411 503
                size = len(self.molecule.xyzs[loc])
412 504
                Size += size
413 505
                apot_added = np.zeros((size, size))
414 506
                bpot_added = np.zeros((size))
415 -
                for xyz, gridxyz, espval, efval  in zip(self.molecule.xyzs[loc:loc+i], self.molecule.gridxyzs[loc:loc+i], self.molecule.espvals[loc:loc+i], self.molecule.efvals[loc:loc+i]):
507 +
                for mol, xyz, gridxyz, espval, efval  in zip(self.molecule.mols[loc+add:loc+i+add], self.molecule.xyzs[loc:loc+i], self.molecule.gridxyzs[loc:loc+i], self.molecule.espvals[loc:loc+i], self.molecule.efvals[loc:loc+i]):
416 508
                    if 'matrices' in self.inp.restraintinfo:
417 509
                        if self.inp.restraintinfo['matrices'] == ['esp', 'ef']:
418 -
419 510
                            apot, bpot = self.EspEfDesignMatrix(xyz, gridxyz, espval, efval)
420 511
                        elif self.inp.restraintinfo['matrices'] == ['esp']:
421 512
                            apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
422 513
                        elif self.inp.restraintinfo['matrices'] == ['ef']:
423 514
                            apot, bpot = self.EfDesignMatrix(xyz, gridxyz, efval)
515 +
                        elif self.inp.restraintinfo['matrices'] == ['obj2']:
516 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj2')
517 +
                        elif self.inp.restraintinfo['matrices'] == ['obj2_2']:
518 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj2_2')
519 +
                        elif self.inp.restraintinfo['matrices'] == ['obj3']:
520 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj3')
521 +
                        elif self.inp.restraintinfo['matrices'] == ['obj3_2']:
522 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj3_2')
523 +
                        elif self.inp.restraintinfo['matrices'] == ['obj4']:
524 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj4')
525 +
                        elif self.inp.restraintinfo['matrices'] == ['obj4_2']:
526 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj4_2')
527 +
                        elif self.inp.restraintinfo['matrices'] == ['obj5']:
528 +
                            apot, bpot = self.NewEspDesignMatrix(xyz, gridxyz, espval, efval, 'obj5')
424 529
                    else:
425 530
                        apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
426 -
                    # Normalization
427 -
                    print(' max(Aij) before Normalization: ', np.amax(apot) )
428 -
                    #scalingFac = np.amax(apot)
429 -
                    scalingFac = (float(self.inp.gridspace)/0.7)**3
430 -
                    print(' scalingFac : ', scalingFac)
431 -
                    apot *= scalingFac
432 -
                    bpot *= scalingFac
433 -
434 531
                    apot_added += apot
435 532
                    bpot_added += bpot
436 533
                apot_added /= i
437 534
                bpot_added /= i
535 +
                if self.normalization == True:
536 +
                    norm = sci.linalg.norm(apot_added, 'fro')
537 +
                    prefac  = 77.87* len(self.molecule.elems[idx]) / (norm * 6.0)
538 +
                    print('Norm:', sci.linalg.norm(apot_added, 'fro'), ', prefac(77.87*Natom/norm/6):', prefac)
539 +
                    apot_added *= prefac
540 +
                    bpot_added *= prefac
541 +
438 542
                apots.append(apot_added)
439 543
                bpots.append(bpot_added)
440 544
                for j in self.molecule.atomids[idx]:
@@ -442,7 +546,9 @@
Loading
442 546
                for j in self.molecule.elems[idx]:
443 547
                    elem_comb.append(j)
444 548
                loc += i
445 -
549 +
        totN = 0
550 +
        for mol in self.molecule.elems:
551 +
            totN += len(mol)
446 552
        apot_comb  = np.zeros((Size, Size))
447 553
        bpot_comb = np.zeros((Size))
448 554
@@ -459,8 +565,7 @@
Loading
459 565
            for i in range(size):
460 566
                bpot_comb[loc+i] = bpot[i]
461 567
            loc += size
462 -
        # print(apots, apot_comb[0], len(apots[0]),'+',len(apots[1]),'=', len(apot_comb))
463 -
        # input()
568 +
464 569
        return apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb
465 570
466 571
    def combine_chargeinfo(self, listofchargeinfo, nAtoms):
@@ -471,16 +576,12 @@
Loading
471 576
        """
472 577
        loc = 0
473 578
        newlistofchargeinfo = []
474 -
        # print(listofchargeinfo)
475 -
        # input()
476 579
        for idx, chargeinfo in enumerate(listofchargeinfo):
477 580
            for indices, resname, charge in chargeinfo:
478 581
                newindices = [i+loc for i in indices]
479 582
                newchargeinfo = [newindices, resname, charge]
480 583
                newlistofchargeinfo.append(newchargeinfo)
481 584
            loc += nAtoms[idx]
482 -
        # print(newlistofchargeinfo)
483 -
        # input()
484 585
        return newlistofchargeinfo
485 586
486 587
    def get_nAtoms(self, nList):
@@ -520,7 +621,7 @@
Loading
520 621
        equivGroup = self.removeSingleElemList(equivGroup)
521 622
        return equivGroup
522 623
523 -
    def force_symmetry(self, apotInp, bpotInp, elem, atomid, equivGroupInp):
624 +
    def force_symmetry(self, apotInp, bpotInp, elem, atomid, listofweights, equivGroupInp):
524 625
        """
525 626
        Force symmetry on matrices using equivGroupInp provided
526 627
@@ -530,17 +631,13 @@
Loading
530 631
531 632
        elem_sym = copy.deepcopy(elem)
532 633
        atomid_sym = copy.deepcopy(atomid)
533 -
634 +
        weights_sym = copy.deepcopy(listofweights)
534 635
        sym = []
535 -
        # print(equivGroupInp, len(apot_sym), apot_sym.shape)
536 -
        # input()
537 636
        for equivlst in equivGroupInp:
538 637
            for idx, i in enumerate(equivlst):
539 638
                if idx == 0:
540 639
                    a = i
541 640
                elif idx > 0:
542 -
                    # print(a,i)
543 -
                    # print(len(apot_sym[a]), len(apot_sym[i]))
544 641
                    apot_sym[a,:] += apot_sym[i,:]
545 642
                    apot_sym[:,a] += apot_sym[:,i]
546 643
                    bpot_sym[a] += bpot_sym[i]
@@ -552,20 +649,18 @@
Loading
552 649
            bpot_sym = np.delete(bpot_sym, i)
553 650
            elem_sym = np.delete(elem_sym, i)
554 651
            atomid_sym = np.delete(atomid_sym, i)  # after forcing symmetry ,it makes singular matrix problem:/
652 +
            weights_sym = np.delete(weights_sym, i)
555 653
654 +
        return apot_sym, bpot_sym, elem_sym, atomid_sym, weights_sym
556 655
557 -
        return apot_sym, bpot_sym, elem_sym, atomid_sym
558 -
559 -
    def apply_set_charge(self, apotInp, bpotInp, atomidInp, atomidinfo, set_charge):
656 +
    def apply_set_charge(self, apotInp, bpotInp, elemInp, atomidInp, atomidinfo, set_charge, weightsInp, tightness):
560 657
        """
561 658
        fix charges listed in set_charge and calculate linear algebra
562 659
563 660
        """
564 -
565 661
        fixedatoms = []
566 662
        setcharges = []
567 663
        for atom in set_charge:
568 -
            #setcharges.append(set_charge[atom]['charge'])
569 664
            val = {'resname': set_charge[atom]['resname'], 'atomname': set_charge[atom]['atomname']}
570 665
            for idx, atomid in enumerate(atomidInp):
571 666
                if val in atomidinfo[atomid]:
@@ -574,6 +669,8 @@
Loading
574 669
                        setcharges.append(set_charge[atom]['charge'])
575 670
        apot_free = copy.deepcopy(apotInp)
576 671
        bpot_free = copy.deepcopy(bpotInp)
672 +
        elem_free = copy.deepcopy(elemInp)
673 +
        weights_free = copy.deepcopy(weightsInp)
577 674
578 675
        for idx, i in enumerate(fixedatoms):
579 676
            bpot_free -=setcharges[idx]*apot_free[:,i]
@@ -590,6 +687,8 @@
Loading
590 687
            apot_free = np.delete(apot_free, i, axis = 0)
591 688
            apot_free = np.delete(apot_free, i, axis = 1)
592 689
            bpot_free = np.delete(bpot_free, i)
690 +
            elem_free = np.delete(elem_free, i)
691 +
            weights_free = np.delete(weights_free, i)
593 692
        if len(apot_free) == 0:
594 693
            qpot_free = []
595 694
        else:
@@ -619,13 +718,9 @@
Loading
619 718
            loc += i
620 719
            qpots.append(qpot)
621 720
        return qpots
622 -
##########################################################################################################
623 -
###                       Charge fit models (Model2, Model3, two-stg-fit)                              ###
624 -
##########################################################################################################
625 721
626 722
    def write_output(self, mol, qpot, moltype='oemol', outfile= 'out.mol2'):
627 723
        if moltype is 'oemol':
628 -
            # set partial charges
629 724
            for idx, atom in enumerate(mol.GetAtoms()):
630 725
                atom.SetPartialCharge(qpot[idx])
631 726
@@ -636,149 +731,15 @@
Loading
636 731
637 732
        elif moltype is 'rdmol':
638 733
            if isinstance(mol, rdchem.Mol):
639 -
                # what I need:
640 -
                print(' Can not generate mol2 file using RDKit yet, will be implemented soon!')
734 +
                print(' Can not generate mol2 file without using Openeye toolkit yet, will be implemented soon!')
641 735
            else:
642 736
                raise RuntimeError(' Failed to identify the RDKit molecule object!')
643 -
644 -
645 737
        else:
646 738
            if isinstance(mol, Molecule):
647 -
                # Check if given molecule object contains atomtype, resid and resname
648 -
                if 'atomtype' not in mol.Data:
649 -
                    mol.atomtype = mol.elem
650 -
                if 'resname' not in mol.Data:
651 -
                    mol.resname = list('MOL' for i in mol.elem)
652 -
                    mol.resid = list(1 for i in mol.elem)
653 -
                charge_type = 'RESPYTE_CHARGE'
654 -
                out = []
655 -
                for I in selection:
656 -
                    out.append('@<TRIPOS>MOLECULE')
657 -
                    out.append('Mol2 generated by resp_optimizer')
658 -
                    out.append('%i %i 0 0 0' % (len(mol.elem), len(mol.bonds)))
659 -
                    out.append('SMALL')
660 -
                    out.append('%s\n' % charge_type)
661 -
                    out.append('@<TRIPOS>ATOM')
662 -
                    for idx, xyz in enumerate(mol.xyzs[I]):
663 -
                        out.append(format_mol2_atom(idx+1, mol.elem[idx], xyz, mol.atomtype[idx], mol.resid[idx], mol.resname[idx], round(qpot[idx],4)))
664 -
                    out.append('@<TRIPOS>BOND')
665 -
                    for idx, i in enumerate(mol.bonds):
666 -
                        out.append(format_mol2_bond(idx+1, i[0]+1, i[1]+1, 1)) # for now, all bonds are assumed to be single bonds
739 +
                print(' Can not generate mol2 file without using Openeye toolkit yet, will be implemented soon!')
667 740
            else:
668 741
                raise RuntimeError('Failed to identify the molecule object.')
669 742
670 -
    def Model2qpot(self, weight):
671 -
        """
672 -
        Solves for q vector in Aq=B (Model 2) with a given restraint weight(a). No restraint on hydrogen atoms.
673 -
674 -
        Parameters
675 -
        ----------
676 -
        weight : float
677 -
            a, scaled facotr which define sthe asymptotic limits of the strength of the restrant.
678 -
679 -
        Returns
680 -
        -------
681 -
        list
682 -
683 -
            list of extended q vectors
684 -
        """
685 -
        # Combine A matrices and B vectors into apot_comb, bpot_comb
686 -
        apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb = self.combineMatrices()
687 -
688 -
        # THen apply Model 2 restraint  -> newapot_comb
689 -
        newapot_comb = copy.deepcopy(np.array(apot_comb))
690 -
        N = len(newapot_comb)
691 -
        for i in range(N):
692 -
            if elem_comb[i] == 1:
693 -
                pass
694 -
            else:
695 -
                newapot_comb[i][i] += weight
696 -
697 -
        # nAtoms : list of number of atoms in each molecule
698 -
        nAtoms = self.get_nAtoms(self.molecule.elems)
699 -
        # newlistofchareginfo : charge info with 'global indices'
700 -
        newlistofchargeinfo = self.combine_chargeinfo(self.molecule.listofchargeinfo, nAtoms)
701 -
702 -
        # Lagrange Charge Constraint on comb matrices
703 -
        apot_constrained, bpot_constrained = self.LagrangeChargeConstraint( newapot_comb, bpot_comb, newlistofchargeinfo)
704 -
705 -
        # Force symmetry based on the atomid
706 -
        equivGroup = self.get_equivGroup(atomid_comb)
707 -
        # print(equivGroup)
708 -
        # input()
709 -
        apot_sym, bpot_sym, elem_sym, atomid_sym = self.force_symmetry(apot_constrained, bpot_constrained, elem_comb, atomid_comb, equivGroup)
710 -
711 -
        # consider set_charge
712 -
        qpot_sym = self.apply_set_charge(apot_sym, bpot_sym, atomid_sym, self.molecule.atomidinfo, self.inp.set_charge)
713 -
        indices_sym = self.getCondensedIndices(len(apot_comb), equivGroup)
714 -
        qpot_expanded = self.Expandqpot(qpot_sym, indices_sym)
715 -
716 -
        # split qpot_Expanded into qpots for each molecule
717 -
        qpots = self.cut_qpot_comb(qpot_expanded, nAtoms)
718 -
719 -
        # write mol2 files with fitted charges.
720 -
        writeMol2 = False
721 -
        path = os.getcwd()
722 -
        writeMol2 = False
723 -
        path = os.getcwd()
724 -
        if os.path.isdir('%s/resp_output' % path):
725 -
            print(' resp_output dir already exists!!! will overwrite anyway:/')
726 -
            writeMol2 = True
727 -
        else:
728 -
            writeMol2 = True
729 -
            os.mkdir('%s/resp_output' % path)
730 -
        # write text file for collecting result:)
731 -
        with open('%s/resp_output/result.txt' % path,'w') as f:
732 -
            f.write('respyte result\n')
733 -
            f.write(' Model 2')
734 -
            if self.inp.restraintinfo['matrices'] == ['esp', 'ef']:
735 -
                f.write(' RESPF\n')
736 -
            elif self.inp.restraintinfo['matrices'] == ['esp']:
737 -
                f.write(' RESP\n')
738 -
            elif self.inp.restraintinfo['matrices'] == ['ef']:
739 -
                f.write(' REF\n')
740 -
            f.write(' weight = %8.4f\n' % (weight))
741 -
            for idx, qpot in enumerate(qpots):
742 -
                f.write('mol%d\n' % (idx+1))
743 -
                for charge in qpot:
744 -
                    f.write('%8.4f\n' % round(charge,4))
745 -
        print()
746 -
        print('-------------------------------------------------------')
747 -
        print()
748 -
        print('                   Model 2 with a =%.4f' % (weight))
749 -
        loc = 0
750 -
        for idx, i in enumerate(self.molecule.nmols):
751 -
            config = 1
752 -
            for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
753 -
                Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
754 -
                apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
755 -
                print()
756 -
                print('           -=#=-  Molecule %d config. %d -=#=-' % (idx+1, config))
757 -
                print()
758 -
                row_format = "{:>10}" * (5)
759 -
                firstrow = ['no.','atomname','resname','atomid','q(opt)']
760 -
                print(row_format.format(*firstrow))
761 -
                for index, q in enumerate(qpots[idx]):
762 -
                    resname = self.molecule.resnames[idx][index]
763 -
                    atomname = self.molecule.atomnames[idx][index]
764 -
                    atomid = self.molecule.atomids[idx][index]
765 -
                    print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
766 -
                print()
767 -
                print(' espRRMS : ', "%.4f" % self.espRRMS(apot, bpot, qpots[idx], espval))
768 -
                print(' efRRMS  : ', "%.4f" % self.efRRMS(Aef, Bef, qpots[idx], efval))
769 -
                if writeMol2 is True:
770 -
                    if self.inp.cheminformatics == 'openeye':
771 -
                        molt = 'oemol'
772 -
                    elif self.inp.cheminformatics == 'rdkit':
773 -
                        molt = 'rdmol'
774 -
                    else:
775 -
                        molt = 'fbmol'
776 -
                    self.write_output(self.molecule.mols[idx], qpots[idx], moltype = molt,
777 -
                                    outfile = '%s/resp_output/mol%d_conf%d.mol2'% (path, idx+1, config))
778 -
                config += 1
779 -
            loc += i
780 -
        return 0
781 -
782 743
    def Model3Amatrix(self, apotInp, weight, tightness, qpotInp, listofelem):
783 744
        """
784 745
        Builds Model 3 A matrix from A0. No restraint on hydrogen atoms.
@@ -805,40 +766,46 @@
Loading
805 766
        newapot = copy.deepcopy(np.array(apotInp))
806 767
        N =  len(newapot)
807 768
        if N != len(listofelem):
808 -
            print(' List of elements should have the same size with A0 matrix.')
769 +
            print('List of elements should have the same size with A0 matrix.')
809 770
            return False
771 +
        list_of_weights = []
772 +
        for i in range(N):
773 +
            list_of_weights.append(weight)
774 +
810 775
        for i in range(N):
811 776
            if listofelem[i] == 'H' or listofelem[i] == 1:
812 777
                continue
813 778
            else:
814 779
                newapot[i][i] += weight / np.sqrt(qpotInp[i]**2 + tightness**2)
815 -
        return newapot
816 -
817 -
    def Model3qpotFn(self, nAtoms, apot_comb, bpot_comb, weight, tightness, elem_comb, atomid_comb, atomidinfo, chargeinfo_comb, set_charge, equivGroupInp):
818 780
781 +
        return newapot, list_of_weights
819 782
783 +
    def Model3qpotFn(self, nAtoms, apot_comb, bpot_comb, weight, tightness, elem_comb, atomid_comb, atomidinfo, chargeinfo_comb, set_charge, equivGroupInp):
820 784
821 785
        def Model3Iteration(qpot_temp):
822 -
            newapot = self.Model3Amatrix(apot_comb, weight, tightness, qpot_temp, elem_comb)
823 -
            # print(len(elem_comb), len(apot_comb), len(newapot))
824 -
            # input()
786 +
            newapot, list_of_weights = self.Model3Amatrix(apot_comb, weight, tightness, qpot_temp, elem_comb)
825 787
            apot_constrained, bpot_constrained = self.LagrangeChargeConstraint(newapot, bpot_comb, chargeinfo_comb)
826 -
            # Force symmetry based on the atomid
827 -
            #equivGroup = self.get_equivGroup(atomid_comb)
828 -
            apot_sym, bpot_sym, elem_sym, atomid_sym = self.force_symmetry(apot_constrained, bpot_constrained, elem_comb, atomid_comb, equivGroupInp)
829 -
            # consider set_charge
830 -
            qpot_sym = self.apply_set_charge(apot_sym, bpot_sym, atomid_sym, atomidinfo, set_charge)
788 +
            apot_sym, bpot_sym, elem_sym, atomid_sym, weights_sym = self.force_symmetry(apot_constrained, bpot_constrained, elem_comb, atomid_comb, list_of_weights, equivGroupInp)
789 +
            qpot_sym = self.apply_set_charge(apot_sym, bpot_sym, elem_sym, atomid_sym, atomidinfo, set_charge, weights_sym, tightness)
831 790
            indices_sym = self.getCondensedIndices(len(apot_comb), equivGroupInp)
832 791
            qpot_expanded = self.Expandqpot(qpot_sym, indices_sym)
833 -
            return qpot_expanded
792 +
            return qpot_expanded, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym
834 793
835 794
        Size = len(apot_comb)
836 795
        qpotInitial = np.zeros((Size))
837 -
        for i in range(10):
838 -
            qpotNext = Model3Iteration(qpotInitial)
796 +
        for i in range(50):
797 +
            qpotNext, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym = Model3Iteration(qpotInitial)
839 798
            if np.linalg.norm(qpotNext-qpotInitial) < 1e-8: break
840 799
            qpotInitial = qpotNext.copy()
841 800
801 +
        #calculate condition number after fitting
802 +
        Options = {'weights':list_of_weights, 'tightness': tightness}
803 +
        cond1 = self.Ncond( apot_comb, elem_comb, qpotNext, Options)
804 +
        cond2 = self.Ncond( apot_sym, elem_sym, qpot_sym, Options)
805 +
        print('###############################################################')
806 +
        print('  After expansion:', cond1)
807 +
        print('  Before expansion:', cond2)
808 +
        print('###############################################################')
842 809
        # split qpot_Expanded into qpots for each molecule
843 810
        qpots = self.cut_qpot_comb(qpotNext, nAtoms)
844 811
        return qpots
@@ -852,7 +819,6 @@
Loading
852 819
                                  elem_comb, atomid_comb, self.molecule.atomidinfo,
853 820
                                  newlistofchargeinfo, self.inp.set_charge, equivGroup)
854 821
855 -
856 822
        # write mol2 files with fitted charges.
857 823
        writeMol2 = False
858 824
        path = os.getcwd()
@@ -882,13 +848,16 @@
Loading
882 848
        print()
883 849
        print('              Model 3 with a = %.4f' % (weight))
884 850
        loc = 0
851 +
        esprrmss = []
852 +
        efrrmss  = []
885 853
        add = 0
886 854
        for idx, i in enumerate(self.molecule.nmols):
887 855
            if i == 0:
888 856
                add += 1
857 +
                esprrmss.append(0)
858 +
                efrrmss.append(0)
889 859
            else:
890 860
                config = 1
891 -
892 861
                for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
893 862
                    Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
894 863
                    apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
@@ -904,8 +873,13 @@
Loading
904 873
                        atomid = self.molecule.atomids[idx][index]
905 874
                        print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
906 875
                    print()
907 -
                    print(' espRRMS : ', "%.4f" % self.espRRMS(apot, bpot, qpots[idx], espval))
908 -
                    print(' efRRMS  : ', "%.4f" % self.efRRMS(Aef, Bef, qpots[idx], efval))
876 +
                    esprrms = self.espRRMS(apot, bpot, qpots[idx], espval)
877 +
                    efrrms  = self.efRRMS(Aef, Bef, qpots[idx], efval)
878 +
                    esprrmss.append(esprrms)
879 +
                    efrrmss.append(efrrms)
880 +
881 +
                    print(' espRRMS : ', "%.4f" % esprrms)
882 +
                    print(' efRRMS  : ', "%.4f" % efrrms)
909 883
                    if writeMol2 is True:
910 884
                        if self.inp.cheminformatics == 'openeye':
911 885
                            molt = 'oemol'
@@ -916,14 +890,13 @@
Loading
916 890
                        self.write_output(self.molecule.mols[loc+config-1+ add], qpots[idx], moltype = molt,                                                                                            outfile = '%s/resp_output/mol%d_conf%d.mol2'% (path, idx+1, config))
917 891
                    config += 1
918 892
            loc += i
919 -
        return 0
893 +
        return qpots, esprrmss, efrrmss
920 894
921 895
    def twoStageFit(self, weight1, weight2, tightness):
922 896
923 -
        apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb = self.combineMatrices()
897 +
        apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb= self.combineMatrices()
924 898
        nAtoms = self.get_nAtoms(self.molecule.elems)
925 -
        equivGroup_comb = self.get_equivGroup(atomid_comb) # equivGroup_combined
926 -
        # I dont think this 'nonpolars' is that useful. But for now, I need it.
899 +
        equivGroup_comb = self.get_equivGroup(atomid_comb)
927 900
        nonpolars = []
928 901
        for index, listofpolar in enumerate(self.molecule.listofpolars):
929 902
            nonpolar = []
@@ -933,7 +906,6 @@
Loading
933 906
                else:
934 907
                    nonpolar.append(i)
935 908
            nonpolars.append(nonpolar)
936 -
        # print(nonpolars)
937 909
        listofpolar_comb = []
938 910
        listofnonpolar_comb = []
939 911
        loc = 0
@@ -952,11 +924,6 @@
Loading
952 924
                listofPolarEquiv.append(equivGroup)
953 925
            else:
954 926
                listofNonpolarEquiv.append(equivGroup)
955 -
        # print('listofpolar_comb', listofpolar_comb)
956 -
        # print('listofnonpolar_comb', listofnonpolar_comb)
957 -
        # print('atomid_comb', atomid_comb)
958 -
        # print(equivGroup_comb)
959 -
        # input()
960 927
        """
961 928
        1st stage: force symmetry on polar atoms and apply smaller restraint weight.
962 929
@@ -965,6 +932,49 @@
Loading
965 932
        qpots_stg1 = self.Model3qpotFn(nAtoms, apot_comb, bpot_comb, weight1, tightness,
966 933
                                       elem_comb, atomid_comb, self.molecule.atomidinfo,
967 934
                                       newlistofchargeinfo, self.inp.set_charge, listofPolarEquiv)
935 +
936 +
        print()
937 +
        print('-------------------------------------------------------')
938 +
        print()
939 +
        print('     Charges from the first stage (for checking)        ')
940 +
        loc = 0
941 +
        esprrmss = []
942 +
        efrrmss  = []
943 +
        add = 0
944 +
        for idx, i in enumerate(self.molecule.nmols):
945 +
            if i == 0:
946 +
                add += 1
947 +
                esprrmss.append(0)
948 +
                efrrmss.append(0)
949 +
            else:
950 +
                config = 1
951 +
                for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
952 +
                    Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
953 +
                    apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
954 +
955 +
                    print()
956 +
                    print('   -=#=- (!!1st stg!!) Molecule %d config. %d -=#=-' % (idx+1, config))
957 +
                    print()
958 +
                    row_format = "{:>10}" * (5)
959 +
                    firstrow = ['no.','atomname','resname','atomid','q(opt)']
960 +
                    print(row_format.format(*firstrow))
961 +
                    for index, q in enumerate(qpots_stg1[idx]):
962 +
                        resname = self.molecule.resnames[idx][index]
963 +
                        atomname = self.molecule.atomnames[idx][index]
964 +
                        atomid = self.molecule.atomids[idx][index]
965 +
                        print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
966 +
                    print()
967 +
                    esprrms = self.espRRMS(apot, bpot, qpots_stg1[idx], espval)
968 +
                    efrrms  = self.efRRMS(Aef, Bef, qpots_stg1[idx], efval)
969 +
                    esprrmss.append(esprrms)
970 +
                    efrrmss.append(efrrms)
971 +
                    print(' espRRMS(1st stg) : ', "%.4f" % esprrms )
972 +
                    print(' efRRMS(1st stg)  : ', "%.4f" % efrrms  )
973 +
                    config += 1
974 +
975 +
            loc += i
976 +
977 +
968 978
        """
969 979
        2nd stage: force symmetry on nonpolar atoms and while fixing fitted polar charge from 1st stg.
970 980
@@ -976,7 +986,6 @@
Loading
976 986
            for charge in qpot:
977 987
                qpot_nonpolar.append(charge)
978 988
979 -
        # cal. nonpolar charge
980 989
        nonpolarchargeinfo = []
981 990
        for idx, chargeinfo in enumerate(newlistofchargeinfo):
982 991
            indices, resname, charge = chargeinfo
@@ -988,23 +997,16 @@
Loading
988 997
                else:
989 998
                    nonpolarindices.append(listofnonpolar_comb.index(i))
990 999
            nonpolarchargeinfo.append([nonpolarindices,resname, charge-polarcharges])
991 -
        # print(newlistofchargeinfo, nonpolarchargeinfo)
992 -
        #input()
993 -
994 1000
        atomid_nonpolar = [atomid_comb[i] for i in listofnonpolar_comb]
995 1001
        elem_nonpolar = [elem_comb[i] for i in listofnonpolar_comb]
996 1002
997 1003
        newlistofNonpolarEquiv  = []
998 -
        # print('#########################################################')
999 -
        # print('listofNonpolarEquiv', listofNonpolarEquiv)
1000 -
        # print('listofnonpolar_comb', listofnonpolar_comb)
1001 1004
        for nonpolarequivgroup in listofNonpolarEquiv:
1002 1005
            newindices = []
1003 1006
            for i in nonpolarequivgroup:
1004 1007
                idx = listofnonpolar_comb.index(i)
1005 1008
                newindices.append(idx)
1006 1009
            newlistofNonpolarEquiv.append(newindices)
1007 -
        # print('newlistofNonpolarEquiv', newlistofNonpolarEquiv)
1008 1010
        listofpolar_comb.sort()
1009 1011
        for i in listofpolar_comb:
1010 1012
            bpot_nonpolar -= qpot_nonpolar[i]*apot_nonpolar[:,i]
@@ -1018,7 +1020,6 @@
Loading
1018 1020
        for idx, Natom in enumerate(nAtoms):
1019 1021
            Nnonpolar = Natom - len(self.molecule.listofpolars[idx])
1020 1022
            Nnonpolars.append(Nnonpolar)
1021 -
        # print('Nnonpolar', Nnonpolar)
1022 1023
        qpots_stg2 = self.Model3qpotFn(Nnonpolars, apot_nonpolar, bpot_nonpolar, weight2, tightness,
1023 1024
                                  elem_nonpolar, atomid_nonpolar, self.molecule.atomidinfo,
1024 1025
                                  nonpolarchargeinfo, self.inp.set_charge, newlistofNonpolarEquiv)
@@ -1058,16 +1059,375 @@
Loading
1058 1059
        print()
1059 1060
        print('    Two stage fit with a1 = %.4f, a2 = %.4f' % (weight1, weight2))
1060 1061
        loc = 0
1062 +
        esprrmss = []
1063 +
        efrrmss  = []
1064 +
        add = 0
1065 +
        for idx, i in enumerate(self.molecule.nmols):
1066 +
            if i == 0:
1067 +
                add += 1
1068 +
                esprrmss.append(0)
1069 +
                efrrmss.append(0)
1070 +
            else:
1071 +
                config = 1
1072 +
                for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
1073 +
                    Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
1074 +
                    apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
1075 +
1076 +
                    print()
1077 +
                    print('           -=#=-  Molecule %d config. %d -=#=-' % (idx+1, config))
1078 +
                    print()
1079 +
                    row_format = "{:>10}" * (5)
1080 +
                    firstrow = ['no.','atomname','resname','atomid','q(opt)']
1081 +
                    print(row_format.format(*firstrow))
1082 +
                    for index, q in enumerate(qpots[idx]):
1083 +
                        resname = self.molecule.resnames[idx][index]
1084 +
                        atomname = self.molecule.atomnames[idx][index]
1085 +
                        atomid = self.molecule.atomids[idx][index]
1086 +
                        print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
1087 +
                    print('MM dipole from two-stg-fit: ', self.MMdipole(xyz, self.molecule.elems[idx], qpots[idx]))
1088 +
                    print()
1089 +
                    esprrms = self.espRRMS(apot, bpot, qpots[idx], espval)
1090 +
                    efrrms  = self.efRRMS(Aef, Bef, qpots[idx], efval)
1091 +
                    esprrmss.append(esprrms)
1092 +
                    efrrmss.append(efrrms)
1093 +
                    print(' espRRMS : ', "%.4f" % esprrms)
1094 +
                    print(' efRRMS  : ', "%.4f" % efrrms)
1095 +
                    if writeMol2 is True:
1096 +
                        if self.inp.cheminformatics == 'openeye':
1097 +
                            molt = 'oemol'
1098 +
                        elif self.inp.cheminformatics == 'rdkit':
1099 +
                            molt = 'rdmol'
1100 +
                        else:
1101 +
                            molt = 'fbmol'
1102 +
                        self.write_output(self.molecule.mols[loc+config-1+ add], qpots[idx], moltype = molt, outfile = '%s/resp_output/mol%d_conf%d.mol2'% (path, idx+1, config))
1103 +
                    config += 1
1104 +
1105 +
            loc += i
1106 +
1107 +
        return  qpots, esprrmss, efrrmss
1108 +
1109 +
    def Model4Amatrix(self, apotInp, weight1, weight2, tightness, qpotInp, listofelem,  listofpolar_comb):
1110 +
        """
1111 +
        Builds Model 4 A matrix from A0. No restraint on hydrogen atoms.
1112 +
1113 +
        Parameters
1114 +
        ----------
1115 +
1116 +
        apotInp : np.ndarray
1117 +
            "A" matrix; 2D array with dimension (# of atoms)
1118 +
        weight1 : float
1119 +
            a1, scaled factor(for nonpolar atoms) which defines the asymptotic limits of the strength of the restraint.
1120 +
        weight2 : float
1121 +
            a2, scaled factor(for polar atoms) which defines the asymptotic limits of the strength of the restraint.
1122 +
        tightness : float
1123 +
            b, tightness of the hyperbola around its minimum.
1124 +
        qpotInp : list
1125 +
            'q' vector
1126 +
        listofelem : list
1127 +
            list of elements whose charges are being fitted
1128 +
1129 +
        Returns
1130 +
        -------
1131 +
        np.ndarray
1132 +
            Model 4 A matrix
1133 +
        """
1134 +
        newapot = copy.deepcopy(np.array(apotInp))
1135 +
        N =  len(newapot)
1136 +
        list_of_weights = []
1137 +
        for i in range(N):
1138 +
            if i in listofpolar_comb:
1139 +
                list_of_weights.append(weight2)
1140 +
            else:
1141 +
                list_of_weights.append(weight1)
1142 +
1143 +
        if N != len(listofelem):
1144 +
            print('List of elements should have the same size with A0 matrix.')
1145 +
            return False
1146 +
        for i in range(N):
1147 +
            if listofelem[i] == 'H' or listofelem[i] == 1:
1148 +
                continue
1149 +
            else:
1150 +
                newapot[i][i] += list_of_weights[i] / np.sqrt(qpotInp[i]**2 + tightness**2)
1151 +
        return newapot, list_of_weights
1152 +
1153 +
    def Model4qpotFn(self, nAtoms, apot_comb, bpot_comb, weight1, weight2, tightness, elem_comb, atomid_comb,
1154 +
                     atomidinfo, chargeinfo_comb, set_charge, equivGroupInp, listofpolar_comb):
1155 +
1156 +
        def Model4Iteration(qpot_temp):
1157 +
            newapot, list_of_weights = self.Model4Amatrix(apot_comb, weight1, weight2, tightness, qpot_temp, elem_comb,  listofpolar_comb)
1158 +
1159 +
            apot_constrained, bpot_constrained = self.LagrangeChargeConstraint(newapot, bpot_comb, chargeinfo_comb)
1160 +
            apot_sym, bpot_sym, elem_sym, atomid_sym, weights_sym= self.force_symmetry(apot_constrained, bpot_constrained, elem_comb, atomid_comb, list_of_weights, equivGroupInp)
1161 +
            # consider set_charge
1162 +
            qpot_sym = self.apply_set_charge(apot_sym, bpot_sym, elem_sym, atomid_sym, atomidinfo, set_charge, weights_sym, tightness)
1163 +
            indices_sym = self.getCondensedIndices(len(apot_comb), equivGroupInp)
1164 +
            qpot_expanded = self.Expandqpot(qpot_sym, indices_sym)
1165 +
            return qpot_expanded, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym
1166 +
1167 +
        Size = len(apot_comb)
1168 +
        qpotInitial = np.zeros((Size))
1169 +
        for i in range(50):
1170 +
            qpotNext, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym = Model4Iteration(qpotInitial)
1171 +
            if np.linalg.norm(qpotNext-qpotInitial) < 1e-8: break
1172 +
            qpotInitial = qpotNext.copy()
1173 +
        #calculate condition number after fitting
1174 +
        Options = {'weights':list_of_weights, 'tightness': tightness}
1175 +
        cond1 = self.Ncond( apot_comb, elem_comb, qpotNext, Options)
1176 +
        cond2 = self.Ncond( apot_sym, elem_sym, qpot_sym, Options)
1177 +
        print('###############################################################')
1178 +
        print('  After expansion:', cond1)
1179 +
        print('  Before expansion:', cond2)
1180 +
        print('###############################################################')
1181 +
        # split qpot_Expanded into qpots for each molecule
1182 +
        qpots = self.cut_qpot_comb(qpotNext, nAtoms)
1183 +
        return qpots
1184 +
1185 +
    def Model4qpot(self, weight1, weight2, tightness):
1186 +
        apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb = self.combineMatrices()
1187 +
        nAtoms = self.get_nAtoms(self.molecule.elems)
1188 +
        equivGroup = self.get_equivGroup(atomid_comb)
1189 +
1190 +
        nonpolars = []
1191 +
        for index, listofpolar in enumerate(self.molecule.listofpolars):
1192 +
            nonpolar = []
1193 +
            for i in range(nAtoms[index]):
1194 +
                if i in listofpolar:
1195 +
                    pass
1196 +
                else:
1197 +
                    nonpolar.append(i)
1198 +
            nonpolars.append(nonpolar)
1199 +
        # print(nonpolars)
1200 +
        listofpolar_comb = []
1201 +
        listofnonpolar_comb = []
1202 +
        loc = 0
1203 +
        for idx, listofpolar in enumerate(self.molecule.listofpolars):
1204 +
            for polaridx in listofpolar:
1205 +
                listofpolar_comb.append(loc+polaridx)
1206 +
            loc += len(self.molecule.elems[idx])
1207 +
        for i in range(len(elem_comb)):
1208 +
            if i not in listofpolar_comb:
1209 +
                listofnonpolar_comb.append(i)
1210 +
1211 +
        newlistofchargeinfo = self.combine_chargeinfo(self.molecule.listofchargeinfo, nAtoms)
1212 +
        qpots = self.Model4qpotFn(nAtoms, apot_comb, bpot_comb, weight1, weight2, tightness,
1213 +
                                  elem_comb, atomid_comb, self.molecule.atomidinfo,
1214 +
                                  newlistofchargeinfo, self.inp.set_charge, equivGroup, listofpolar_comb)
1215 +
1216 +
1217 +
        # write mol2 files with fitted charges.
1218 +
        writeMol2 = False
1219 +
        path = os.getcwd()
1220 +
        if os.path.isdir('%s/resp_output' % path):
1221 +
            print(' resp_output dir already exists!!! will overwrite anyway:/')
1222 +
            writeMol2 = True
1223 +
        else:
1224 +
            writeMol2 = True
1225 +
            os.mkdir('%s/resp_output' % path)
1226 +
        # write text file for collecting result:)
1227 +
        with open('%s/resp_output/result.txt' % path,'w') as f:
1228 +
            f.write('respyte result\n')
1229 +
            f.write(' Model 4')
1230 +
            if self.inp.restraintinfo['matrices'] == ['esp', 'ef']:
1231 +
                f.write(' RESPF\n')
1232 +
            elif self.inp.restraintinfo['matrices'] == ['esp']:
1233 +
                f.write(' RESP\n')
1234 +
            elif self.inp.restraintinfo['matrices'] == ['ef']:
1235 +
                f.write(' REF\n')
1236 +
            f.write(' weight1 = %8.4f, weight2 = %8.4f, tightness = %8.4f\n' % (weight1, weight2, tightness))
1237 +
            for idx, qpot in enumerate(qpots):
1238 +
                f.write('mol%d\n' % (idx+1))
1239 +
                for charge in qpot:
1240 +
                    f.write('%8.4f\n' % round(charge,4))
1241 +
        print()
1242 +
        print('-------------------------------------------------------')
1243 +
        print()
1244 +
        print('              Model 4 with a1 = %.4f, a2 = %.4f' % (weight1, weight2))
1245 +
        loc = 0
1246 +
        esprrmss = []
1247 +
        efrrmss  = []
1061 1248
        add = 0
1062 1249
        for idx, i in enumerate(self.molecule.nmols):
1063 1250
            if i == 0:
1064 1251
                add += 1
1252 +
                esprrmss.append(0)
1253 +
                efrrmss.append(0)
1065 1254
            else:
1066 1255
                config = 1
1067 1256
                for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
1068 1257
                    Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
1069 1258
                    apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
1259 +
                    print()
1260 +
                    print('           -=#=-  Molecule %d config. %d -=#=-' % (idx+1, config))
1261 +
                    print()
1262 +
                    row_format = "{:>10}" * (5)
1263 +
                    firstrow = ['no.','atomname','resname','atomid','q(opt)']
1264 +
                    print(row_format.format(*firstrow))
1265 +
                    for index, q in enumerate(qpots[idx]):
1266 +
                        resname = self.molecule.resnames[idx][index]
1267 +
                        atomname = self.molecule.atomnames[idx][index]
1268 +
                        atomid = self.molecule.atomids[idx][index]
1269 +
                        print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
1270 +
                    print()
1271 +
                    print('MM dipole from model4: ', self.MMdipole(xyz, self.molecule.elems[idx], qpots[idx]))
1272 +
                    print()
1273 +
                    esprrms = self.espRRMS(apot, bpot, qpots[idx], espval)
1274 +
                    efrrms  = self.efRRMS(Aef, Bef, qpots[idx], efval)
1275 +
                    esprrmss.append(esprrms)
1276 +
                    efrrmss.append(efrrms)
1277 +
                    print(' espRRMS : ', "%.4f" % esprrms)
1278 +
                    print(' efRRMS  : ', "%.4f" % efrrms)
1279 +
                    if writeMol2 is True:
1280 +
                        if self.inp.cheminformatics == 'openeye':
1281 +
                            molt = 'oemol'
1282 +
                        elif self.inp.cheminformatics == 'rdkit':
1283 +
                            molt = 'rdmol'
1284 +
                        else:
1285 +
                            molt = 'fbmol'
1286 +
                        self.write_output(self.molecule.mols[loc+config-1+ add], qpots[idx], moltype = molt,                                                                                            outfile = '%s/resp_output/mol%d_conf%d.mol2'% (path, idx+1, config))
1287 +
                    config += 1
1288 +
            loc += i
1289 +
        return qpots, esprrmss, efrrmss
1290 +
1291 +
    def Model6Amatrix(self, apotInp, weight, tightness, qpotInp, listofelem,  listofburied):
1292 +
        """
1293 +
        Builds Model 4 A matrix from A0. No restraint on hydrogen atoms.
1294 +
1295 +
        Parameters
1296 +
        ----------
1297 +
1298 +
        apotInp : np.ndarray
1299 +
            "A" matrix; 2D array with dimension (# of atoms)
1300 +
        weight1 : float
1301 +
            a1, scaled factor(for buried atoms) which defines the asymptotic limits of the strength of the restraint.
1302 +
        weight2 : float
1303 +
            a2, scaled factor(for exposed atoms) which defines the asymptotic limits of the strength of the restraint.
1304 +
        tightness : float
1305 +
            b, tightness of the hyperbola around its minimum.
1306 +
        qpotInp : list
1307 +
            'q' vector
1308 +
        listofelem : list
1309 +
            list of elements whose charges are being fitted
1310 +
1311 +
        Returns
1312 +
        -------
1313 +
        np.ndarray
1314 +
            Model 6 A matrix
1315 +
        """
1316 +
        newapot = copy.deepcopy(np.array(apotInp))
1317 +
        N =  len(newapot)
1318 +
        if N != len(listofelem):
1319 +
            print('List of elements should have the same size with A0 matrix.')
1320 +
            return False
1321 +
        
1322 +
        list_of_weights = []
1323 +
        for i in range(N):
1324 +
            if i in listofburied:
1325 +
                list_of_weights.append(weight)
1326 +
            else:
1327 +
                list_of_weights.append(weight/10)
1328 +
        for i in range(N):
1329 +
            if listofelem[i] == 'H' or listofelem[i] == 1:
1330 +
                continue
1331 +
            else:
1332 +
                newapot[i][i] += list_of_weights[i]/ np.sqrt(qpotInp[i]**2 + tightness**2)
1333 +
        return newapot, list_of_weights
1334 +
1335 +
    def Model6qpotFn(self, nAtoms, apot_comb, bpot_comb, weight, tightness, elem_comb, atomid_comb,
1336 +
                     atomidinfo, chargeinfo_comb, set_charge, equivGroupInp,  listofburied_comb):
1337 +
1338 +
        def Model6Iteration(qpot_temp):
1339 +
            newapot, list_of_weights = self.Model6Amatrix(apot_comb, weight, tightness, qpot_temp, elem_comb,   listofburied_comb)
1340 +
            apot_constrained, bpot_constrained = self.LagrangeChargeConstraint(newapot, bpot_comb, chargeinfo_comb)
1341 +
            # Force symmetry based on the atomid
1342 +
            apot_sym, bpot_sym, elem_sym, atomid_sym, weights_sym = self.force_symmetry(apot_constrained, bpot_constrained, elem_comb, atomid_comb, list_of_weights, equivGroupInp)
1343 +
            # consider set_charge
1344 +
            qpot_sym= self.apply_set_charge(apot_sym, bpot_sym, elem_sym, atomid_sym, atomidinfo, set_charge, weights_sym, tightness)
1345 +
            indices_sym = self.getCondensedIndices(len(apot_comb), equivGroupInp)
1346 +
            qpot_expanded = self.Expandqpot(qpot_sym, indices_sym)
1347 +
            return qpot_expanded, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym
1348 +
1349 +
        Size = len(apot_comb)
1350 +
        qpotInitial = np.zeros((Size))
1351 +
        for i in range(50):
1352 +
            qpotNext, list_of_weights, apot_sym, bpot_sym, elem_sym, qpot_sym = Model6Iteration(qpotInitial)
1353 +
            if np.linalg.norm(qpotNext-qpotInitial) < 1e-8: break
1354 +
            qpotInitial = qpotNext.copy()
1355 +
1356 +
        #calculate condition number after fitting
1357 +
        Options = {'weights':list_of_weights, 'tightness': tightness}
1358 +
        cond1 = self.Ncond( apot_comb, elem_comb, qpotNext, Options)
1359 +
        cond2 = self.Ncond( apot_sym, elem_sym, qpot_sym, Options)
1360 +
        print('###############################################################')
1361 +
        print('  After expansion:', cond1)
1362 +
        print('  Before expansion:', cond2)
1363 +
        print('###############################################################')
1364 +
        # split qpot_Expanded into qpots for each molecule
1365 +
        qpots = self.cut_qpot_comb(qpotNext, nAtoms)
1366 +
        return qpots
1070 1367
1368 +
    def Model6qpot(self, weight, tightness):
1369 +
        apots, bpots, atomid_comb, elem_comb, apot_comb, bpot_comb = self.combineMatrices()
1370 +
        nAtoms = self.get_nAtoms(self.molecule.elems)
1371 +
        equivGroup = self.get_equivGroup(atomid_comb)
1372 +
1373 +
        listofburied_comb = []
1374 +
        loc = 0
1375 +
        for idx, elems in enumerate(self.molecule.elems):
1376 +
            if len(elems) == 0:
1377 +
                continue
1378 +
            else:
1379 +
                for atmidx in self.molecule.listofburieds[idx]:
1380 +
                    listofburied_comb.append(atmidx+loc)
1381 +
                loc += len(elems)
1382 +
1383 +
        newlistofchargeinfo = self.combine_chargeinfo(self.molecule.listofchargeinfo, nAtoms)
1384 +
        qpots = self.Model6qpotFn(nAtoms, apot_comb, bpot_comb, weight, tightness,
1385 +
                                  elem_comb, atomid_comb, self.molecule.atomidinfo,
1386 +
                                  newlistofchargeinfo, self.inp.set_charge, equivGroup, listofburied_comb)
1387 +
1388 +
1389 +
        # write mol2 files with fitted charges.
1390 +
        writeMol2 = False
1391 +
        path = os.getcwd()
1392 +
        if os.path.isdir('%s/resp_output' % path):
1393 +
            print(' resp_output dir already exists!!! will overwrite anyway:/')
1394 +
            writeMol2 = True
1395 +
        else:
1396 +
            writeMol2 = True
1397 +
            os.mkdir('%s/resp_output' % path)
1398 +
        # write text file for collecting result:)
1399 +
        with open('%s/resp_output/result.txt' % path,'w') as f:
1400 +
            f.write('respyte result\n')
1401 +
            f.write(' Model 6')
1402 +
            if self.inp.restraintinfo['matrices'] == ['esp', 'ef']:
1403 +
                f.write(' RESPF\n')
1404 +
            elif self.inp.restraintinfo['matrices'] == ['esp']:
1405 +
                f.write(' RESP\n')
1406 +
            elif self.inp.restraintinfo['matrices'] == ['ef']:
1407 +
                f.write(' REF\n')
1408 +
            f.write(' weight = %8.4f, tightness = %8.4f\n' % (weight, tightness))
1409 +
            for idx, qpot in enumerate(qpots):
1410 +
                f.write('mol%d\n' % (idx+1))
1411 +
                for charge in qpot:
1412 +
                    f.write('%8.4f\n' % round(charge,4))
1413 +
        print()
1414 +
        print('-------------------------------------------------------')
1415 +
        print()
1416 +
        print('              Model 6 with a = %.4f' % (weight))
1417 +
        loc = 0
1418 +
        esprrmss = []
1419 +
        efrrmss  = []
1420 +
        add = 0
1421 +
        for idx, i in enumerate(self.molecule.nmols):
1422 +
            if i == 0:
1423 +
                add += 1
1424 +
                esprrmss.append(0)
1425 +
                efrrmss.append(0)
1426 +
            else:
1427 +
                config = 1
1428 +
                for xyz,gridxyz, espval, efval  in zip(self.molecule.xyzs[loc: loc+i], self.molecule.gridxyzs[loc: loc+i], self.molecule.espvals[loc: loc+i], self.molecule.efvals[loc: loc+i]):
1429 +
                    Aef, Bef = self.EfDesignMatrix(xyz, gridxyz, efval)
1430 +
                    apot, bpot = self.EspDesignMatrix(xyz, gridxyz, espval)
1071 1431
                    print()
1072 1432
                    print('           -=#=-  Molecule %d config. %d -=#=-' % (idx+1, config))
1073 1433
                    print()
@@ -1080,8 +1440,14 @@
Loading
1080 1440
                        atomid = self.molecule.atomids[idx][index]
1081 1441
                        print(row_format.format(index, atomname, resname, atomid,"%.4f" % q))
1082 1442
                    print()
1083 -
                    print(' espRRMS : ', "%.4f" % self.espRRMS(apot, bpot, qpots[idx], espval))
1084 -
                    print(' efRRMS  : ', "%.4f" % self.efRRMS(Aef, Bef, qpots[idx], efval))
1443 +
                    print('MM dipole from model6: ', self.MMdipole(xyz, self.molecule.elems[idx], qpots[idx]))
1444 +
                    print()
1445 +
                    esprrms = self.espRRMS(apot, bpot, qpots[idx], espval)
1446 +
                    efrrms  = self.efRRMS(Aef, Bef, qpots[idx], efval)
1447 +
                    esprrmss.append(esprrms)
1448 +
                    efrrmss.append(efrrms)
1449 +
                    print(' espRRMS : ', "%.4f" % esprrms)
1450 +
                    print(' efRRMS  : ', "%.4f" % efrrms)
1085 1451
                    if writeMol2 is True:
1086 1452
                        if self.inp.cheminformatics == 'openeye':
1087 1453
                            molt = 'oemol'
@@ -1092,4 +1458,51 @@
Loading
1092 1458
                        self.write_output(self.molecule.mols[loc+config-1+ add], qpots[idx], moltype = molt,                                                                                            outfile = '%s/resp_output/mol%d_conf%d.mol2'% (path, idx+1, config))
1093 1459
                    config += 1
1094 1460
            loc += i
1095 -
        return 0
1461 +
        return qpots, esprrmss, efrrmss
1462 +
1463 +
    def EspResiduals(self, qpotsInp):
1464 +
        list_of_residuals = []
1465 +
        if len(molecule.xyzs) == len(qpotsInp):
1466 +
            qpots = qpotsInp
1467 +
        else:
1468 +
            qpots = []
1469 +
            for idxx, Nconf in enumerate(molecule.nmols):
1470 +
                for k in range(Nconf):
1471 +
                    qpots.append(qpotsInp[idxx])
1472 +
            if len(molecule.xyzs) != len(qpots):
1473 +
                raise RuntimeError('sth wrong!!')
1474 +
        for idx, xyzs in enumerate(self.molecule.xyzs):
1475 +
            molbohr = xyzs
1476 +
            gridxyzs = self.molecule.gridxyzs[idx]
1477 +
            espval = self.molecule.espvals[idx]
1478 +
1479 +
            x = np.array([xyz[0] for xyz in molbohr])
1480 +
            y = np.array([xyz[1] for xyz in molbohr])
1481 +
            z = np.array([xyz[2] for xyz in molbohr])
1482 +
            potv = np.array(espval)
1483 +
            potvsq = potv*potv
1484 +
            ssvpot = np.sum(potvsq)
1485 +
            sspotsorted = np.sum(np.sort(potvsq))
1486 +
            invRij = []
1487 +
            dxij = []
1488 +
            dyij = []
1489 +
            dzij = []
1490 +
            for pt in range(len(gridxyzs)):
1491 +
                i = pt
1492 +
                dxj = x - gridxyzs[i,0]
1493 +
                dxij.append(dxj)
1494 +
                dyj = y - gridxyzs[i,1]
1495 +
                dyij.append(dyj)
1496 +
                dzj = z - gridxyzs[i,2]
1497 +
                dzij.append(dzj)
1498 +
                rijsq = dxj*dxj + dyj*dyj + dzj*dzj
1499 +
                invRjSq = 1.0/rijsq
1500 +
                invRij.append(np.sqrt(invRjSq))
1501 +
1502 +
            vi = np.zeros(len(gridxyzs))
1503 +
            for j, invRj in enumerate(invRij):
1504 +
                vj = qpots[idx]*invRj
1505 +
                vi[j] = np.sum(np.sort(vj))
1506 +
            residual = vi - espval
1507 +
            list_of_residuals.append(residual)
1508 +
        return list_of_residuals

@@ -4,7 +4,7 @@
Loading
4 4
# from molecule_resp import Molecule_HJ
5 5
from collections import OrderedDict
6 6
from warnings import warn
7 -
from .molecule import *
7 +
from respyte.molecule import *
8 8
9 9
class Engine(object):
10 10
    def __init__(self, coordfile = None, input_file = None):

@@ -1,29 +1,34 @@
Loading
1 1
from respyte.molecule import *
2 2
3 3
BondiRadii = [1.2, 1.4, # exchanged None to 2.0
4 -
              1.81, 2.0, 2.0, 1.70, 1.55, 1.52, 1.47, 1.54,
5 -
              2.27, 1.73, 2.0, 2.22, 1.80, 1.80, 1.75, 1.88,
6 -
              2.75, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.63, 1.40, 1.39, 1.87, 2.0, 1.85, 1.90, 1.83, 2.02]
4 +
              1.81, 2.00, 2.00, 1.70, 1.55, 1.52, 1.47, 1.54,
5 +
              2.27, 1.73, 2.00, 2.22, 1.80, 1.80, 1.75, 1.88,
6 +
              2.75, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 1.63, 1.40, 1.39, 1.87, 2.00, 1.85, 1.90, 1.83, 2.02,
7 +
              2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 1.63, 1.72, 1.62, 1.93, 2.17, 2.00, 2.00, 1.98, 2.16]
7 8
# from oechem
8 -
BondiRadiiOechem = [1.2, 1.4,
9 -
                    1.82, 2.0, 2.0, 1.70, 1.55, 1.52, 1.47, 1.54,
10 -
                    2.27, 1.73, 2.0, 2.1, 1.80, 1.80, 1.75, 1.88,
11 -
                    2.75, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.63, 1.40, 1.39, 1.87, 2.0, 1.85, 1.90, 1.85, 2.02]
9 +
BondiRadiiOechem = [1.20, 1.40,
10 +
                    1.82, 2.00, 2.00, 1.70, 1.55, 1.52, 1.47, 1.54,
11 +
                    2.27, 1.73, 2.00, 2.01, 1.80, 1.80, 1.75, 1.88,
12 +
                    2.75, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 1.63, 1.40, 1.39, 1.87, 2.00, 1.85, 1.90, 1.85, 2.02,
13 +
                    2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 1.72, 1.58, 2.00, 2.17, 2.00, 2.06, 1.98, 2.16]
12 14
13 15
AlvarezRadii = [1.20, 1.43,
14 16
                2.12, 1.98, 1.91, 1.77, 1.66, 1.50, 1.46, 1.58,
15 17
                2.50, 2.51, 2.25, 2.19, 1.90, 1.89, 1.82, 1.83,
16 -
                2.73, 2.62, 2.58, 2.46, 2.42, 2.45, 2.45, 2.44, 2.40, 2.40, 2.38, 2.39, 2.32, 2.29, 1.88, 1.82, 1.86, 2.25]
18 +
                2.73, 2.62, 2.58, 2.46, 2.42, 2.45, 2.45, 2.44, 2.40, 2.40, 2.38, 2.39, 2.32, 2.29, 1.88, 1.82, 1.86, 2.25,
19 +
                3.21, 2.84, 2.75, 2.52, 2.56, 2.45, 2.44, 2.46, 2.44, 2.15, 2.53, 2.49, 2.43, 2.42, 2.47, 1.99, 2.04, 2.06]
17 20
PolarHradii = 0.95 # in Angstrom
18 21
bohr2Ang = 0.52918825
19 22
20 23
21 24
def SelectGridPts(mol, inner, outer, pts, radiiType):
22 25
    if radiiType == 'bondi':
23 -
        radiiType = BondiRadii
26 +
        # radiiType = BondiRadii
27 +
        radiiType = BondiRadiiOechem
24 28
    elif radiiType == 'Alvarez':
25 29
        radiiType = AlvarezRadii
26 30
    culled = []
31 +
    selectedPts = []
27 32
    xyzs = mol.xyzs[0]
28 33
    innersSq = []
29 34
    outersSq = []
@@ -33,9 +38,6 @@
Loading
33 38
        radii = float(radiiType[int(elem)-1])
34 39
        innersSq.append(radii*radii*inner*inner)
35 40
        outersSq.append(radii*radii*outer*outer)
36 -
    # print('inners:', innersSq)
37 -
    # print('outers: ', outersSq)
38 -
    # input()
39 41
    for idx, pt in enumerate(pts):
40 42
        goodPt = False
41 43
        for xyz, innerSq, outerSq in zip(xyzs, innersSq, outersSq):
@@ -50,5 +52,7 @@
Loading
50 52
                break
51 53
        if goodPt:
52 54
            culled.append(idx)
55 +
            selectedPts.append(pt)
56 +
57 +
    return culled, selectedPts
53 58
54 -
    return culled

@@ -34,6 +34,8 @@
Loading
34 34
            self.restraintinfo = {}
35 35
            self.gridinfo = None
36 36
            self.gridspace = 0.7
37 +
            self.normalization = False
38 +
            self.symmetry  = True
37 39
        else:
38 40
            self.readinp(inputFile)
39 41
@@ -76,21 +78,13 @@
Loading
76 78
            newcharge_equal.append([atomname, resname])
77 79
        charge_equal = newcharge_equal
78 80
79 -
        # Check how many molecules and how many conformers are set to be fitted
80 81
        nmols = []
81 82
        for mol in inp['molecules']:
82 83
            nmols.append(int(inp['molecules'][mol]))
83 84
84 -
        # Read charge fit model (Model2, Model3, 2-stg-fit)
85 85
        restraintinfo = inp['restraint']
86 86
87 87
        gridinfo = {}
88 -
        # if 'grid_info' in inp:
89 -
        #     gridinfo['type'] = inp['grid_info']['type']
90 -
        #     gridinfo['radii'] = inp ['grid_info']['radii']
91 -
92 -
        # if 'shell_select' in inp:
93 -
        #    gridinfo['shell_select'] = inp['shell_select']
94 88
        if 'boundary_select' in inp:
95 89
            gridinfo['boundary_select'] = inp['boundary_select']
96 90
@@ -99,6 +93,16 @@
Loading
99 93
        else:
100 94
            space = 0.7
101 95
96 +
        if 'normalization' in inp:
97 +
            normalization = inp['normalization']
98 +
        else:
99 +
            normalization = False
100 +
101 +
        if 'symmetry' in inp:
102 +
            symmetry = inp['symmetry']
103 +
        else:
104 +
            symmetry = True
105 +
102 106
        self.cheminformatics = cheminformatics
103 107
        self.set_charge      = set_charge
104 108
        self.resChargeDict   = resChargeDict
@@ -107,10 +111,12 @@
Loading
107 111
        self.restraintinfo   = restraintinfo
108 112
        self.gridinfo        = gridinfo
109 113
        self.gridspace       = space
110 -
114 +
        self.normalization   = normalization
115 +
        self.symmetry        = symmetry
111 116
def main():
112 117
    inp = Input()
113 -
    inp.readinput(sys.argv[1])
118 +
    inp.readinp(sys.argv[1])
119 +
    print(inp.normalization)
114 120
115 121
if __name__ == "__main__":
116 122
    main()
Files Coverage
respyte 8.68%
Project Totals (14 files) 8.68%
6.1
TRAVIS_OS_NAME=osx
6.3
TRAVIS_PYTHON_VERSION=3.5
TRAVIS_OS_NAME=linux
6.4
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
6.2
TRAVIS_OS_NAME=osx
1
# Codecov configuration to make it a bit less noisy
2
coverage:
3
  status:
4
    patch: false
5
    project:
6
      default:
7
        threshold: 50%
8
comment:
9
  layout: "header"
10
  require_changes: false
11
  branches: null
12
  behavior: default
13
  flags: null
14
  paths: null
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading