| 1 | # -*- coding: utf-8 -*-
|
|---|
| 2 | '''
|
|---|
| 3 | Created on Aug 31, 2010
|
|---|
| 4 |
|
|---|
| 5 | @author: Jan H. Meinke
|
|---|
| 6 | '''
|
|---|
| 7 |
|
|---|
| 8 | import smmp_p as smmp
|
|---|
| 9 | import numpy as np
|
|---|
| 10 | from math import sqrt, cos
|
|---|
| 11 |
|
|---|
| 12 | def scan():
|
|---|
| 13 | """Iterate through all interaction pairs and return lists of
|
|---|
| 14 | interaction partners for each type of interaction.
|
|---|
| 15 |
|
|---|
| 16 | This implementation is almost equivalent to the Fortran
|
|---|
| 17 | implementation in enyshe.f.
|
|---|
| 18 |
|
|---|
| 19 | """
|
|---|
| 20 | cutoff = 32
|
|---|
| 21 | longDistancePartners = {}
|
|---|
| 22 | oneFourPartners = {}
|
|---|
| 23 |
|
|---|
| 24 | EvdW = 0
|
|---|
| 25 | Ec = 0
|
|---|
| 26 | Ehb = 0
|
|---|
| 27 | Etor = 0
|
|---|
| 28 | ifivr = smmp.mol_i.ivrml1[0]
|
|---|
| 29 | ntlvr = smmp.mol_par.nvr
|
|---|
| 30 | i1s = smmp.mol_i.imsml1[0] + smmp.mol_i.nmsml[0]
|
|---|
| 31 | # Loop over variables in reverse order.
|
|---|
| 32 | for iv in smmp.var_i.iorvr[ifivr + ntlvr - 2:ifivr:-1] - 1:
|
|---|
| 33 | ia = smmp.var_i.iatvr[iv] - 1
|
|---|
| 34 | ic = smmp.var_i.iclvr[iv] - 1
|
|---|
| 35 | if smmp.var_i.ityvr[iv] == 3:
|
|---|
| 36 | Etor += smmp.epar_r.e0to[ic] * (1.0 + smmp.epar_r.sgto[ic] *
|
|---|
| 37 | cos(smmp.atm_r.toat[ia] * smmp.epar_r.rnto[ic]))
|
|---|
| 38 |
|
|---|
| 39 | i2s = i1s - 1
|
|---|
| 40 | i1s = smmp.var_i.imsvr1[iv - 1]
|
|---|
| 41 | for i1, i2 in zip(smmp.var_i.latms1[i1s:i2s] - 1,
|
|---|
| 42 | smmp.var_i.latms2[i1s:i2s]):
|
|---|
| 43 | for i, ity, cqi, xi, yi, zi, ivw1, ivw2, i141, i142 in zip(
|
|---|
| 44 | np.arange(i1, i2), smmp.atm_i.ityat[i1:i2],
|
|---|
| 45 | smmp.epar_r.conv * smmp.atm_r.cgat[i1:i2],
|
|---|
| 46 | smmp.atm_r.xat[i1:i2],
|
|---|
| 47 | smmp.atm_r.yat[i1:i2], smmp.atm_r.zat[i1:i2],
|
|---|
| 48 | smmp.eny_i.ivwat1[i1:i2] - 1, smmp.eny_i.ivwat2[i1:i2],
|
|---|
| 49 | smmp.eny_i.i14at1[i1:i2] - 1,smmp.eny_i.i14at1[i1:i2]):
|
|---|
| 50 | longDistancePartners[i] = []
|
|---|
| 51 | oneFourPartners[i] = []
|
|---|
| 52 | for j, jty, xij, yij, zij, cqj in zip(np.arange(ivw1, ivw2),
|
|---|
| 53 | smmp.atm_i.ityat[ivw1:ivw2],
|
|---|
| 54 | smmp.atm_r.xat[ivw1:ivw2] - xi,
|
|---|
| 55 | smmp.atm_r.yat[ivw1:ivw2] - yi,
|
|---|
| 56 | smmp.atm_r.zat[ivw1:ivw2] - zi,
|
|---|
| 57 | smmp.atm_r.cgat[ivw1:ivw2]):
|
|---|
| 58 | if abs(i - j) < cutoff:
|
|---|
| 59 | longDistancePartners[i].append(j)
|
|---|
| 60 | rij2 = xij * xij + yij * yij + zij * zij
|
|---|
| 61 | rij4 = rij2 * rij2
|
|---|
| 62 | rij6 = rij4 * rij2
|
|---|
| 63 | rij = sqrt(rij2)
|
|---|
| 64 | Ec += cqi * cqj / rij
|
|---|
| 65 | if smmp.epar_i.ihbty[ity, jty]:
|
|---|
| 66 | Ehb += smmp.epar_r.ahb[ity, jty] / (rij6 * rij6) -\
|
|---|
| 67 | smmp.epar_r.chb[ity, jty] / (rij6 * rij4)
|
|---|
| 68 | else:
|
|---|
| 69 | EvdW += smmp.epar_r.aij[ity, jty] / (rij6 * rij6) -\
|
|---|
| 70 | smmp.epar_r.cij[ity, jty] / rij
|
|---|
| 71 |
|
|---|
| 72 | for j in smmp.eny_i.l14at[i141:i142]:
|
|---|
| 73 | jty = smmp.atm_i.ityat[j]
|
|---|
| 74 | xij = smmp.atm_r.xat[j] - xi
|
|---|
| 75 | yij = smmp.atm_r.xat[j] - yi
|
|---|
| 76 | zij = smmp.atm_r.xat[j] - zi
|
|---|
| 77 | rij2 = xij * xij + yij * yij + zij * zij
|
|---|
| 78 | rij4 = rij2 * rij2
|
|---|
| 79 | rij6 = rij4 * rij2
|
|---|
| 80 | rij = sqrt(rij2)
|
|---|
| 81 | Ec += cqi * cqj / rij
|
|---|
| 82 | if smmp.epar_i.ihbty[ity, jty]:
|
|---|
| 83 | Ehb += smmp.epar_r.ahb[ity, jty] / (rij6 * rij6) -\
|
|---|
| 84 | smmp.epar_r.chb[ity, jty] / (rij6 * rij4)
|
|---|
| 85 | longDistancePartners[i].append(j)
|
|---|
| 86 | else:
|
|---|
| 87 | oneFourPartners[i].append(j)
|
|---|
| 88 | EvdW += smmp.epar_r.a14[ity, jty] / (rij6 * rij6) -\
|
|---|
| 89 | smmp.epar_r.cij[ity, jty] / rij
|
|---|
| 90 |
|
|---|
| 91 | if not longDistancePartners[i]:
|
|---|
| 92 | del(longDistancePartners[i])
|
|---|
| 93 | if not oneFourPartners[i]:
|
|---|
| 94 | del(oneFourPartners[i])
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | print EvdW + Ec + Ehb + Etor, EvdW, Ec, Ehb, Etor
|
|---|
| 98 |
|
|---|
| 99 | return longDistancePartners, oneFourPartners
|
|---|
| 100 |
|
|---|
