[SciPy-User] Molecular Index Calculation

[Hjalmar Turesson]

Hi Stephen,


I see that N = 10 (ie. the length of Z), and that r is indexed by i and j,
which run up to N. However, r is only a 5x5 array. Is this correct?

Otherwise, something like this might work:

import numpy as np

start=1
finish=31
points=300
s = np.linspace(start, finish, points)

MASS =  np.array([12.011, 1.008, 79.9, 35.453, 18.998])

r = np.array([[0., 2.059801, 3.60937686, 3.32591826, 2.81569212],
              [2.059801, 0., 4.71452879, 4.45776445, 4.00467382],
              [3.60937686, 4.71452879, 0., 5.66500917, 5.26602175],
              [3.32591826, 4.45776445, 5.66500917, 0., 5.02324896],
              [2.81569212, 4.00467382, 5.26602175, 5.02324896, 0.]])

r[r == 0.] = 0.0001

Z = np.array([6.0, 210.0, 102.0, 54.0, 35.0, 17.0, 9.0, 595.0, 315.0,
153.0])

N = Z.size

I = np.zeros(s.shape)

for i in range(1, N):
    for j in range(j):
        I += Z[i] * Z[j] * np.sin(s * r[i, j])/(s * r[i, j])


Do you have an example of the correct output?

Best,
Hjalmar

On Sat, Apr 15, 2017 at 6:29 AM, Stephen P. Molnar <s.molnar at sbcglobal.net>
wrote:

> I am a newcomer to Python and am attempting to translate a FORTRAN program
> that I wrote about 20 years ago into Python, specifically v-3.6.
>
> Let me say that I am not asking someone to write the program for me, but
> only to point me is the right direction.
>
> So far, I have bumbled my way to the point that I can get all of the input
> data resulting from a quantum mechanical calculation of a very simple
> organic molecule in to a Python program, but am encountering problems with
> processing the data.
>
> The equation that I a want to evaluate (the one that I programmed in
> FORTRAN) is equation (7) in the attached scan of one of the pages of the
> following document.
>
> Stephen P. Molnar and James W. King, Theory and Applications of the
> Integrated Molecular Transform and the Normalized Molecular Moment
> Structure Descriptors: QSAR and QSPR Paradigms, Int. J Quantum Chem., 85,
> 662 (2001).
>
> the equation is
>
> I(s) = sum from i=2toN sum from j=1toi-1 Z_sub_i*Z_sub_j *
>
>    sin  s*r_sub_ij
>   -----------------  Equation (7)
>        s*r_sub_ij
>
> What I have managed to do so far is for a simple organic molecule
> containing 5 atoms (CHFClBr):
>
> import numpy as np
> import numpy, pandas, scipy.spatial.distance as dist
>
> r=[]        #distances between pairs of atoms
> z = []      #atomic numbers of atoms
> Z_dist = [] #products of successive atomic numbers
> I = []      #vector calculated in equation  (7)
>
> start=1
> finish=31
> points=300
> s = np.linspace(start,finish,points)
>
>
> name = input("Enter Molecule ID: ")
> name = str(name)
> name_in = name+'.dat'
>
> df = pandas.read_table(name_in, skiprows=2, sep=" ", skipinitialspace=True)
> z = numpy.array(df["ZA"])
> N =  numpy.ma.size(z)  ## of atoms in molecule
> a = numpy.array(df[["X", "Y", "Z"]])
> dist.squareform(dist.pdist(a, "euclidean"))
> anrows, ancols = np.shape(a)
> a_new = a.reshape(anrows, 1, ancols)
> diff = a_new - a
>
> r = (diff ** 2).sum(2)
> r = np.sqrt(r)
>
> for j in range(0,N-1):
>      for k in range(j+1,N):
>          Z_diff = (z[j]*z[k])
>
>
> For the test molecule I get the following:
>
> MASS =  [ 12.011   1.008  79.9    35.453  18.998]  (Atomic weights of
> atoms)
>
> r:  [ 0.          2.059801    3.60937686  3.32591826  2.81569212]
>      [ 2.059801    0.          4.71452879  4.45776445  4.00467382]
>      [ 3.60937686  4.71452879  0.          5.66500917  5.26602175]
>      [ 3.32591826  4.45776445  5.66500917  0.          5.02324896]
>      [ 2.81569212  4.00467382  5.26602175  5.02324896 0.        ]
>
> z:
>
> 6.0
> 210.0
> 102.0
> 54.0
> 35.0
> 17.0
> 9.0
> 595.0
> 315.0
> 153.0
>
> I have checked these calculations with a spreadsheet calculation,
> consequently I'm correct as far as I have gotten.
>
> However, my problem is calculating the value of I for the 300 x-coordinate
> values required for the molecular index calculation.
>
> The problem that I am currently having is calculating the product of the
> pre-sine function and the sine term for the 300 values of s.
>
> Unfortunately, the wealth of  function in Python makes it difficult for me
> to ascertain just how to proceed.  It seems that I have the problem of
> seeing the trees because of the size of the forest.  At this point Goggling
> only makes my confusion deeper.  Any pointers in the correct direction will
> be greatly appreciated.
>
> Thanks in advance.
>
> --
> Stephen P. Molnar, Ph.D.                Life is a fuzzy set
> www.molecular-modeling.net              Stochastic and multivariate
> (614)312-7528 (c)
> Skype: smolnar1
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at python.org
> https://mail.python.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-user/attachments/20170415/d8987020/attachment.html>