#!/usr/bin/python3.6

from numpy import *

import string
import random
import math
import readandwrite
import intcoord
import takeinp
import bmatrix
import mymath
import datastruct
import various
import physconstants
import constraints
import int_values
import sys

'''
program vibFreq.py
	- determines hessian and dynamical matrix and writes their eigenvalues and eigenfunctions 
	in several file formats

	- files *.molden can be visualized using molden, jmol, or other visualization software
	supporting molden frequency and normal coordinates format

	- file hesseModes.dat contains eigenvalues and eigenvectors of Hesse matrix sorted in 
	ascending order by the eigenvalues

	- example of usage:
 
	python3 ../../tools/vibFreq/vibFreq.py 0.2 0.5 0.1 2 2 2 > freq.dat

        in this case, fractional coordinates are shifted by 0.2 0.5 and 0.1 and the cell is 
	multiplied by 2 along each lattice axis. Motions corresponding to different frequencies 
	and force constants obtained for such a supercell are written in files 
	vibModes_shifted.molden and hesseModes_shifted.molden, respectively (the analogous files 
	without "shifted" in their name contain a similar information for unshifted unit cell).
'''

def find_translations(lattmat, rmax):
    lattinv = linalg.inv(transpose(lattmat))
    translations = [1, 1, 1]
    for i in range(3):
        translations[i] = int(round(rmax*sum(lattinv[i]**2)**0.5))+1
    return translations


def dynMat_(h, m, ftags):
    mmat = zeros((len(h), len(h)), float)
    indx = 0
    for i in range(len(m)):
        if ftags[3*i] == 0:
            mmat[indx][indx] = m[i]**(-0.5)
            indx += 1
        if ftags[3*i+1] == 0:
            mmat[indx][indx] = m[i]**(-0.5)
            indx += 1
        if ftags[3*i+2] == 0:
            mmat[indx][indx] = m[i]**(-0.5)
            indx += 1
    dynmat = dot(mmat, h)
    dynmat = dot(dynmat, mmat)
    return dynmat


def dynMat(h, m):
    mmat = zeros((len(h), len(h)), float)
    for i in range(len(m)):
        mmat[3*i][3*i] = m[i]**(-0.5)
        mmat[3*i+1][3*i+1] = m[i]**(-0.5)
        mmat[3*i+2][3*i+2] = m[i]**(-0.5)
    dynmat = dot(mmat, h)
    dynmat = dot(dynmat, mmat)
    return dynmat


def read_forces(file):
    f = open(file, 'r')
    index = 0
    gradients = []
    positions = []
    step = 0.
    for line in f.readlines():
        if index > 0:
            index += 1
        prefield = line.split()
        if len(prefield) > 0:
            if prefield[0] == 'POTIM':
                step = float(prefield[2])
            if prefield[0] == 'POSITION':
                cartgrad = []
                cartpos = []
                index = 1
            if index > 2:
                if len(prefield) < 6:
                    gradients.append(array(cartgrad))
                    positions.append(array(cartpos))
                    index = 0
                    continue
                for j in range(0, 3):
                    test = float(prefield[j])
                    cartpos.append(test)
                for j in range(3, 6):
                    test = float(prefield[j])
                    cartgrad.append(test)
    f.close()
    return array(gradients), array(positions)


shift=array([0.5,0.5,0.5])
translations=array([1,1,1])

if len(sys.argv)==4:
    shift[0] = float(sys.argv[1])
    shift[1] = float(sys.argv[2])
    shift[2] = float(sys.argv[3])

if len(sys.argv)==7:
    shift[0] = float(sys.argv[1])
    shift[1] = float(sys.argv[2])
    shift[2] = float(sys.argv[3])
    translations[0] = float(sys.argv[4])
    translations[1] = float(sys.argv[5])
    translations[2] = float(sys.argv[6])

totalTrans=translations[0]*translations[1]*translations[2]


pC = physconstants.physConstants()


amutokg = 1.66053886e-27
evtoj = 1.60217646e-19
angsttom = 1e-10
cl = 2.997925e08
planck = 6.62608e-34
potim = 1.0
au2a = 0.5291772083
temperature = 100
colprob = 0.0

# mass in kg
# mass=readMasses('MASSES')*amutokg
# flags=readFlags('FLAGS')
# print "mass",mass
# print "flags",flags
tknp = takeinp.TakeInput()
tknp.read("OUTCAR")
# tknp.convert_to_au()

numofatoms = tknp.numofatoms
nfree=tknp.nfree
lattmat = tknp.lattmat
katoms = tknp.types
coords_d = tknp.coords_d
coords_c = tknp.coords_c
atomquality = tknp.atomicFlags
attags = []
#print("katoms", katoms)
#print("atomquality", atomquality)
for i in range(len(katoms)):
    attags += katoms[i]*[atomquality[i]]


rmax = 30.
trans = find_translations(lattmat, rmax)
cnexp = [9., 14.]


try:
    f = open('MASSES', 'r')
    mass = []
    for line in f.readlines():
        line = string.split(line)
        if len(line) > 0:
            mass.append(float(line[0]))
    f.close()
    mass = array(mass)
    mass *= amutokg

except IOError:
    mass_ = tknp.atomicMass
    #print(mass_)
    mass = []
    for i in range(len(katoms)):
        mass += katoms[i]*[mass_[i]]
    mass = array(mass)
    mass *= amutokg


allowedAtoms = []
fixedAtoms = []
allowedDirections = ones(3*numofatoms, float)

try:
    f = open('ALLOWED', 'r')
except IOError:
    #print("file with allowed atoms not detected")
    pass
else:
    allowedDirections = zeros(3*numofatoms, float)
    allowedAtoms = []
    for line in f.readlines():
        line = string.split(line)
        for i in range(len(line)):
            dummy = int(line[i])
            if dummy > 0 and dummy <= numofatoms:
                dummy -= 1
                allowedAtoms.append(dummy)
                allowedDirections[3*dummy] = 1.
                allowedDirections[3*dummy+1] = 1.
                allowedDirections[3*dummy+2] = 1.
    f.close()
    allowedAtoms = array(allowedAtoms)
    print(allowedAtoms)

gradients, positions = read_forces('OUTCAR')

# c remove contribution from atoms fixed by user
# for i in range(len(positions)):
#  for j in range(3*numofatoms):
#    positions[i,j]*=allowedDirections[j]

if nfree==2:
    bmat = zeros((int((len(gradients)-1)/2), 3*numofatoms),float)
    hes1 = zeros((int((len(gradients)-1)/2), 3*numofatoms),float)
    for i in range(int((len(gradients)-1)/2)):
        stepDirection = (positions[2*i+1]-positions[2*i+2])*allowedDirections
        # step=sum((positions[2*i+1]-positions[2*i+2])**2)**0.5
        step = (sum(stepDirection*stepDirection))**0.5
        # print i,step
        if step > 1e-5:
            bmat[i] = (positions[2*i+1]-positions[2*i+2])/step  # (2*step)
            hes1[i] = (gradients[2*i+2]-gradients[2*i+1])/step  # (2*step)
        else:
            bmat[i] = 0.
            hes1[i] = 0.
elif nfree==1:
    bmat = zeros(((len(gradients)-1), 3*numofatoms), float)
    hes1 = zeros(((len(gradients)-1), 3*numofatoms), float)
    for i in range(1, len(gradients)):
        stepDirection = (positions[i]-positions[0])*allowedDirections
        step = (sum(stepDirection*stepDirection))**0.5
        if step > 1e-5:
            bmat[i-1] = (positions[i]-positions[0])/step  # (2*step)
            hes1[i-1] = -(gradients[i]-gradients[0])/step  # (2*step)
        else:
            bmat[i] = 0.
            hes1[i] = 0.
else:
    sys.exit("unsupported NFREE detected: "+str(nfree))   

hessian = dot(transpose(bmat), hes1)

# c determine the fixed atoms
fcount = 0
tcount = 0
ftags = zeros(3*numofatoms, int)
for i in range(len(hessian)):
    if sum(abs(hessian[i])) < 1e-3:
        ftags[i] = 1
        fcount += 1
    else:
        tcount += 1

# if fcount>0:
#  hessian_=[]
#  for i in range(len(hessian)):
#    if ftags[i]==0:
#      hline=[]
#      for j in range(len(hessian)):
#        if ftags[j]==0:
#          hline.append(hessian[i][j])
#      hessian_.append(hline)
#  hessian=array(hessian_)

# print len(hessian),len(hessian[0])
# print dagf

# fcount=0;tcount=0
# for i in range(len(hessian)):
#  print i,sum(abs(hessian[i]))
#  if sum(abs(hessian[i]))<1e-3:
#    fcount+=1
#  else:
#    tcount+=1
#
# print fcount, tcount
# print das

hessian = (hessian+transpose(hessian))/2

# c take care of fixed atoms
# for i in range(len(hessian)):
#  if ftags[i]==1:
#    hessian[i][:]=0.
#    hessian[:][i]=0.

# convert into internals and back to cartesians in order
# to get rid of translations and rotations
# atomquality=read_atomicFlags("OUTCAR")
# ASCALE=1.1;BSCALE=1.1;FRAG=2;RELAX=0;TORSIONS=1;SUBST=100
ATRADII = []
COVALENTRADII = datastruct.Elementprop().covalentradii
for i in range(len(atomquality)):
    index = atomquality[i]
    ATRADII.append(COVALENTRADII[index])

ASCALE = 2.10*pC.AU2A
BSCALE = 1.30
FRAG = 2
RELAX = 0
TORSIONS = 1
SUBST = 100
# ASCALE=1.1*pC.AU2A
intrn = intcoord.Intern(atomquality, ATRADII, ASCALE, BSCALE, FRAG, RELAX,
                        TORSIONS, SUBST, numofatoms, lattmat, coords_d, coords_c, katoms)
# intrn=intcoord.Intern(atomquality,ASCALE,BSCALE,FRAG,RELAX,TORSIONS,SUBST,'POSCAR')
primcoords = intrn.internalcoords
cartesian = various.change_format(coords_c)
#coords_d_shifted = coords_d+array([0.5, 0.5, 0.5])
coords_d_shifted = coords_d+shift
for i in range(numofatoms):
    coords_d_shifted[i, 0] = coords_d_shifted[i, 0] % 1
    coords_d_shifted[i, 1] = coords_d_shifted[i, 1] % 1
    coords_d_shifted[i, 2] = coords_d_shifted[i, 2] % 1

coords_c_shifted = dot(coords_d_shifted, lattmat)

primcoords, CONSTRAINTS, conststat, ircdata = constraints.get_iconst(
    primcoords, cartesian, lattmat, attags)


deal = int_values.Int_values(
    coords_c, primcoords, lattmat, ircdata, cnexp, trans)

b = bmatrix.bmatrix(cartesian, primcoords, numofatoms,
                    lattmat, RELAX, cnexp, ircdata, trans)
Bmat = b.Bmatrix

if len(CONSTRAINTS) > 0:
    Bmat = constraints.setup_iconst(CONSTRAINTS, Bmat)

    for i in CONSTRAINTS:
        Bmat[i] = 0.

# c complement
# for i in range(len(Bmat)):
#  lzero=0
#  for j in CONSTRAINTS:
#    if i==j:lzero=1
#  if lzero==0: Bmat[i]=0.


Gmatrix = dot(transpose(Bmat), Bmat)
GmatrixInv = mymath.mygeneralized_inverse(Gmatrix)
Amat = dot(GmatrixInv, transpose(Bmat))


transmat = mymath.cd_transmatrix(lattmat, len(cartesian))

# if fcount>0:
#  transmat_=[]
#  for i in range(len(transmat)):
#    if ftags[i]==0:
#      hline=[]
#      for j in range(len(transmat)):
#        if ftags[j]==0:
#          hline.append(hessian[i][j])
#      transmat_.append(hline)
#  transmat=array(transmat_)


hessian2 = dot(transpose(transmat), hessian)
hessian2 = dot(hessian2, transmat)
hessian2 = dot(transpose(Amat), hessian2)
hessian2 = dot(hessian2, Amat)

# eigval,eigvect=mymath.eigh_tb(hessian2)
# print eigval

hessian2 = dot(hessian2, (Bmat))
hessian2 = dot(transpose(Bmat), hessian2)
transmat = linalg.inv(transmat)
hessian2 = dot(transpose(transmat), hessian2)
hessian2 = dot(hessian2, transmat)
hessian2 = (hessian2+transpose(hessian2))/2


# c take care of fixed atoms
if fcount > 0:
    hessian_ = []
    for i in range(len(hessian2)):
        if ftags[i] == 0:
            hline = []
            for j in range(len(hessian2)):
                if ftags[j] == 0:
                    hline.append(hessian2[i][j])
            hessian_.append(hline)
    hessian2 = array(hessian_)


eigval, eigvect = mymath.eigh_tb(hessian2)

# c write down the eigevalues and eigenmodes for the Hesse matrix:
f = open('hesseModes.dat', 'w')
for i in range(len(eigval)):
    counth = 0
    rrr = str(i+1)+' f  ='+'\t'+str(eigval[i]) + ' eV/A^2'+'\n'
    f.write(rrr)
    for j in range(numofatoms):
        indx = 3*j
        if ftags[indx] == 1:
            xx = round(0.0, 6)
        else:
            xx = round(eigvect[i][counth], 6)
            counth += 1

        indx = 3*j+1
        if ftags[indx] == 1:
            yy = round(0.0, 6)
        else:
            yy = round(eigvect[i][counth], 6)
            counth += 1

        indx = 3*j+2
        if ftags[indx] == 1:
            zz = round(0.0, 6)
        else:
            zz = round(eigvect[i][counth], 6)
            counth += 1

        rrr = " "+'%.6f' % (xx)
        rrr += " "+'%.6f' % (yy)
        rrr += " "+'%.6f' % (zz)+'\n'
        f.write(rrr)
    f.write('\n')

f.close()

##########################################################################
# c write eigen-values and vectors of Hesse matrix in molden format
coords_c /= au2a
coords_c_shifted /= au2a

f = open('hesseModes.molden', 'w')
f2 = open('hesseModes_shifted.molden', 'w')

rrr = " [Molden Format]"+'\n'+" [FREQ]"+'\n'
f.write(rrr)
f2.write(rrr)

for i in range(len(eigval)):
    rrr = '%.2f' % round(eigval[i], 2)+'\n'
    f.write(rrr)
    f2.write(rrr)

rrr = " [FR-COORD]"+'\n'
f.write(rrr)
f2.write(rrr)

for ii in range(translations[0]):
    for jj in range(translations[1]):
        for kk in range(translations[2]):
            for i in range(numofatoms):
                rrr= attags[i]
                rrr2= attags[i]
                dummy1=coords_c[i] + ii*lattmat[0]/au2a + jj*lattmat[1]/au2a + kk*lattmat[2]/au2a
                dummy2=coords_c_shifted[i] + ii*lattmat[0]/au2a + jj*lattmat[1]/au2a + kk*lattmat[2]/au2a
                rrr+=" "+'%.6f'%(round(dummy1[0],6))
                rrr+=" "+'%.6f'%(round(dummy1[1],6))
                rrr+=" "+'%.6f'%(round(dummy1[2],6))+'\n'
                rrr2+=" "+'%.6f'%(round(dummy2[0],6))
                rrr2+=" "+'%.6f'%(round(dummy2[1],6))
                rrr2+=" "+'%.6f'%(round(dummy2[2],6))+'\n'
                f.write(rrr)
                f2.write(rrr2)


rrr = " [FR-NORM-COORD]"+'\n'
f.write(rrr)
f2.write(rrr)
for i in range(len(eigval)):
    rrr = " vibration "+str(i+1)+'\n'
    f.write(rrr)
    f2.write(rrr)

    for jj in range(totalTrans):
        counth=0
        for j in range(numofatoms):
            indx=3*j
            if ftags[indx]==1:
                xx=round(0.0,6)
            else:
                xx=round(eigvect[i][counth],6)
                counth+=1

            indx=3*j+1
            if ftags[indx]==1:
                yy=round(0.0,6)
            else:
                yy=round(eigvect[i][counth],6)
                counth+=1

            indx=3*j+2
            if ftags[indx]==1:
                zz=round(0.0,6)
            else:
                zz=round(eigvect[i][counth],6)
                counth+=1

            rrr=" "+'%.6f'%(xx)
            rrr+=" "+'%.6f'%(yy)
            rrr+=" "+'%.6f'%(zz)+'\n'
            f.write(rrr)
            f2.write(rrr)


f2.close()
f.close()

# c transform back to A
# eigvect*=au2a
coords_c *= au2a
coords_c_shifted *= au2a

######################################################################









dynmat = dynMat_(hessian2, mass, ftags)
# dynmat=dynMat(hessian2,mass)
dynmat = dynmat*evtoj/angsttom**2


# for i in range(len(dynmat)):
#  if ftags[i]==1:
#    dynmat[i][:]=0.
#    dynmat[:][i]=0.

eigval, eigvect = mymath.eigh_tb(dynmat)
eigval = sort(eigval)
for i in range(len(eigval)):
    if eigval[i] > 0:
        print(i+1, 'f', '=', '\t', (eigval[i]**0.5)/cl/(2*math.pi)*1e-02, 'cm-1', '\t',\
            (eigval[i]**0.5)/1e12, '2PiTHz',\
            '\t', eigval[i]**0.5/(1e12*2*pi), 'THz', '\t', \
            1e03*(eigval[i]**0.5)*planck/evtoj/(2*pi), 'meV')
    else:
        print(i+1, 'f/i', '=', '\t', (abs(eigval[i])**0.5)/cl/(2*math.pi)*1e-02, 'cm-1', '\t',\
            (abs(eigval[i])**0.5)/1e12, '2PiTHz',\
            '\t', abs(eigval[i])**0.5/(1e12*2*pi), 'THz', '\t', \
            1e03*(abs(eigval[i])**0.5)*planck/evtoj/(2*pi), 'meV')

######################################################################################
f = open('vibModes.dat', 'w')
for i in range(len(eigval)):
    counth = 0
    if eigval[i] > 0:
        rrr = str(i+1)+' f  ='+'\t'+str((eigval[i]**0.5)/cl/(2*math.pi)*1e-02) + ' cm-1'+'\t'\
            + str((eigval[i]**0.5)/1e12)+' 2PiTHz'\
            + '\t'+str(eigval[i]**0.5/(1e12*2*pi)) + ' THz' + '\t' \
            + str(1e03*(eigval[i]**0.5)*planck/evtoj/(2*pi))+' meV'+'\n'
        f.write(rrr)
        for j in range(numofatoms):
            indx = 3*j
            if ftags[indx] == 1:
                xx = round(0.0, 6)
            else:
                xx = round(eigvect[i][counth], 6)
                counth += 1

            indx = 3*j+1
            if ftags[indx] == 1:
                yy = round(0.0, 6)
            else:
                yy = round(eigvect[i][counth], 6)
                counth += 1

            indx = 3*j+2
            if ftags[indx] == 1:
                zz = round(0.0, 6)
            else:
                zz = round(eigvect[i][counth], 6)
                counth += 1

            rrr = '\t'+attags[j]
            rrr += " " + str(j+1)
            rrr += " "+'%.6f' % (xx)
            rrr += " "+'%.6f' % (yy)
            rrr += " "+'%.6f' % (zz)+'\n'
            f.write(rrr)
        f.write('\n')
    else:
        eigval[i] = abs(eigval[i])
        rrr = str(i+1)+' f/i ='+'\t'+str((eigval[i]**0.5)/cl/(2*math.pi)*1e-02) + ' cm-1'+'\t'\
            + str((eigval[i]**0.5)/1e12)+' 2PiTHz'\
            + '\t'+str(eigval[i]**0.5/(1e12*2*pi)) + ' THz' + '\t' \
            + str(1e03*(eigval[i]**0.5)*planck/evtoj/(2*pi))+' meV'+'\n'
        f.write(rrr)
        for j in range(numofatoms):
            indx = 3*j
            if ftags[indx] == 1:
                xx = round(0.0, 6)
            else:
                xx = round(eigvect[i][counth], 6)
                counth += 1

            indx = 3*j+1
            if ftags[indx] == 1:
                yy = round(0.0, 6)
            else:
                yy = round(eigvect[i][counth], 6)
                counth += 1

            indx = 3*j+2
            if ftags[indx] == 1:
                zz = round(0.0, 6)
            else:
                zz = round(eigvect[i][counth], 6)
                counth += 1

            rrr = '\t'+attags[j]
            rrr += " " + str(j+1)
            rrr += " "+'%.6f' % (xx)
            rrr += " "+'%.6f' % (yy)
            rrr += " "+'%.6f' % (zz)+'\n'
            f.write(rrr)
        f.write('\n')
f.close()

######################################################################################
# as vibFreq.dat but eigenvectors are mass-weighted, NOT  normalized afterwards!
######################################################################################
# c rescale eigenvectors of dynamical matrix by 1/sqrt(m)
#c and re-normalize
eigvect_ = 1*eigvect
mmat = zeros((len(dynmat), len(dynmat)), float)

indx = 0
for i in range(len(mass)):
    if ftags[3*i] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1
    if ftags[3*i+1] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1
    if ftags[3*i+2] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1

# c masses in amu
mmat *= amutokg**0.5

# for i in range(len(mass)):
#  mmat[3*i][3*i]=mass[i]**(-0.5)
#  mmat[3*i+1][3*i+1]=mass[i]**(-0.5)
#  mmat[3*i+2][3*i+2]=mass[i]**(-0.5)

eigvect = dot(eigvect, mmat)

# c DO NOT normalize after reweigthing!
# for i in range(len(eigvect)):
#  eigvect[i]/=(sum(eigvect[i]**2)**0.5)

f = open('vibModes_MW.dat', 'w')
for i in range(len(eigval)):
    counth = 0
    if eigval[i] > 0:
        rrr = str(i+1)+' f  ='+'\t'+str((eigval[i]**0.5)/cl/(2*math.pi)*1e-02) + ' cm-1'+'\t'\
            + str((eigval[i]**0.5)/1e12)+' 2PiTHz'\
            + '\t'+str(eigval[i]**0.5/(1e12*2*pi)) + ' THz' + '\t' \
            + str(1e03*(eigval[i]**0.5)*planck/evtoj/(2*pi))+' meV'+'\n'
        f.write(rrr)
        for j in range(numofatoms):
            indx = 3*j
            if ftags[indx] == 1:
                xx = round(0.0, 8)
            else:
                xx = round(eigvect[i][counth], 6)
                counth += 1

            indx = 3*j+1
            if ftags[indx] == 1:
                yy = round(0.0, 8)
            else:
                yy = round(eigvect[i][counth], 8)
                counth += 1

            indx = 3*j+2
            if ftags[indx] == 1:
                zz = round(0.0, 8)
            else:
                zz = round(eigvect[i][counth], 8)
                counth += 1

            rrr = '\t'+attags[j]
            rrr += " " + str(j+1)
            rrr += " "+'%.8f' % (xx)
            rrr += " "+'%.8f' % (yy)
            rrr += " "+'%.8f' % (zz)+'\n'
            f.write(rrr)
        f.write('\n')
    else:
        eigval[i] = abs(eigval[i])
        rrr = str(i+1)+' f/i ='+'\t'+str((eigval[i]**0.5)/cl/(2*math.pi)*1e-02) + ' cm-1'+'\t'\
            + str((eigval[i]**0.5)/1e12)+' 2PiTHz'\
            + '\t'+str(eigval[i]**0.5/(1e12*2*pi)) + ' THz' + '\t' \
            + str(1e03*(eigval[i]**0.5)*planck/evtoj/(2*pi))+' meV'+'\n'
        f.write(rrr)
        for j in range(numofatoms):
            indx = 3*j
            if ftags[indx] == 1:
                xx = round(0.0, 8)
            else:
                xx = round(eigvect[i][counth], 8)
                counth += 1

            indx = 3*j+1
            if ftags[indx] == 1:
                yy = round(0.0, 8)
            else:
                yy = round(eigvect[i][counth], 8)
                counth += 1

            indx = 3*j+2
            if ftags[indx] == 1:
                zz = round(0.0, 6)
            else:
                zz = round(eigvect[i][counth], 8)
                counth += 1

            rrr = '\t'+attags[j]
            rrr += " " + str(j+1)
            rrr += " "+'%.8f' % (xx)
            rrr += " "+'%.8f' % (yy)
            rrr += " "+'%.8f' % (zz)+'\n'
            f.write(rrr)
        f.write('\n')
f.close()

eigvect = 1*eigvect_

######################################################################################
#######                      Molden format                                     #######
######################################################################################

# c rescale eigenvectors of dynamical matrix by 1/sqrt(m)
#c and re-normalize
eigvect_ = 1*eigvect
mmat = zeros((len(dynmat), len(dynmat)), float)

indx = 0
for i in range(len(mass)):
    if ftags[3*i] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1
    if ftags[3*i+1] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1
    if ftags[3*i+2] == 0:
        mmat[indx][indx] = mass[i]**(-0.5)
        indx += 1

# for i in range(len(mass)):
#  mmat[3*i][3*i]=mass[i]**(-0.5)
#  mmat[3*i+1][3*i+1]=mass[i]**(-0.5)
#  mmat[3*i+2][3*i+2]=mass[i]**(-0.5)

eigvect = dot(eigvect, mmat)

for i in range(len(eigvect)):
    eigvect[i] /= (sum(eigvect[i]**2)**0.5)

# c write the mass-weighted vibrational modes
# c write down the eigevalues and eigenmodes for the Hesse matrix:
f = open('weightedModes.dat', 'w')
for i in range(len(eigval)):
    counth = 0
    #rrr=str(i+1)+' f  ='+'\t'+str(eigval[i]) +'\n'
    rrr = str(i+1)+' f  ='+'\t' + \
        '%.2f' % round((eigval[i]**0.5)/cl/(2*math.pi)*1e-02, 2)+'\n'
    f.write(rrr)
    for j in range(numofatoms):
        indx = 3*j
        if ftags[indx] == 1:
            xx = round(0.0, 6)
        else:
            xx = round(eigvect[i][counth], 6)
            counth += 1

        indx = 3*j+1
        if ftags[indx] == 1:
            yy = round(0.0, 6)
        else:
            yy = round(eigvect[i][counth], 6)
            counth += 1

        indx = 3*j+2
        if ftags[indx] == 1:
            zz = round(0.0, 6)
        else:
            zz = round(eigvect[i][counth], 6)
            counth += 1

        rrr = " "+'%.6f' % (xx)
        rrr += " "+'%.6f' % (yy)
        rrr += " "+'%.6f' % (zz)+'\n'
        f.write(rrr)
    f.write('\n')

f.close()


# c first transform to atomic unit
# eigvect/=au2a
coords_c /= au2a
coords_c_shifted /= au2a

f = open('vibModes.molden', 'w')
f2 = open('vibModes_shifted.molden', 'w')

rrr = " [Molden Format]"+'\n'+" [FREQ]"+'\n'
f.write(rrr)
f2.write(rrr)
for i in range(len(eigval)):
    if eigval[i] > 0:
        rrr = '%.2f' % round((eigval[i]**0.5)/cl/(2*math.pi)*1e-02, 2)+'\n'
        f.write(rrr)
        f2.write(rrr)
    else:
        rrr = '%.2f' % round(-(abs(eigval[i])
                               ** 0.5)/cl/(2*math.pi)*1e-02, 2)+'\n'
        f.write(rrr)
        f2.write(rrr)
rrr = " [FR-COORD]"+'\n'
f.write(rrr)
f2.write(rrr)

for ii in range(translations[0]):
    for jj in range(translations[1]):
        for kk in range(translations[2]):
            for i in range(numofatoms):
                rrr= attags[i]
                rrr2= attags[i]
                dummy1=coords_c[i] + ii*lattmat[0]/au2a + jj*lattmat[1]/au2a + kk*lattmat[2]/au2a
                dummy2=coords_c_shifted[i] + ii*lattmat[0]/au2a + jj*lattmat[1]/au2a + kk*lattmat[2]/au2a
                rrr+=" "+'%.6f'%(round(dummy1[0],6))
                rrr+=" "+'%.6f'%(round(dummy1[1],6))
                rrr+=" "+'%.6f'%(round(dummy1[2],6))+'\n'
                rrr2+=" "+'%.6f'%(round(dummy2[0],6))
                rrr2+=" "+'%.6f'%(round(dummy2[1],6))
                rrr2+=" "+'%.6f'%(round(dummy2[2],6))+'\n'
                f.write(rrr)
                f2.write(rrr2)

rrr = " [FR-NORM-COORD]"+'\n'
f.write(rrr)
f2.write(rrr)
for i in range(len(eigval)):
    rrr = " vibration "+str(i+1)+'\n'
    f.write(rrr)
    f2.write(rrr)

    for jj in range(totalTrans):
        counth=0
        for j in range(numofatoms):
            indx=3*j
            if ftags[indx]==1:
                xx=round(0.0,6)
            else:
                xx=round(eigvect[i][counth],6)
                counth+=1

            indx=3*j+1
            if ftags[indx]==1:
                yy=round(0.0,6)
            else:
                yy=round(eigvect[i][counth],6)
                counth+=1

            indx=3*j+2
            if ftags[indx]==1:
                zz=round(0.0,6)
            else:
                zz=round(eigvect[i][counth],6)
                counth+=1

            rrr=" "+'%.6f'%(xx)
            rrr+=" "+'%.6f'%(yy)
            rrr+=" "+'%.6f'%(zz)+'\n'
            f.write(rrr)
            f2.write(rrr)

f2.close()
f.close()

# c transform back to A
# eigvect*=au2a
coords_c *= au2a

eigvect = 1*eigvect_

######################################################################################
#######                      a-Climax format                                   #######
######################################################################################

# position in periodic table of elements
pte = {'H': 1, 'He': 2, 'Li': 3, 'Be': 4, 'B': 5, 'C': 6, 'N': 7, 'O': 8, 'F': 9, 'Ne': 10,
       'Na': 11, 'Mg': 12, 'Al': 13, 'Si': 14, 'P': 15, 'S': 16, 'Cl': 17, 'Ar': 18,
       'Cu': 29, 'Ga': 31, 'Br': 35, 'Pd': 46, 'I': 53}

sigmas = {'H': 80., 'He': 0., 'Li': 0.9, 'Be': 2e-3, 'B': 1.7, 'C': 1e-3, 'N': 5e-1, 'O': 8e-3, 'F': 1e-3, 'Ne': 0.,
  'Na': 1.6, 'Mg': 3.7, 'Al': 1e-2, 'Si': 4e-3, 'P': 5e-3, 'S': 1e-2, 'Cl': 5.3, 'Ar': 0,
          'Cu': 0.1, 'Ga': 0.2, 'Br': 0.1, 'Pd': 0.2, 'I': 0.1, }

# c rescale eigenvectors of dynamical matrix by 1/sqrt(m)
#c and re-normalize
#eigvect_ = 1*eigvect
# mmat=zeros((3*numofatoms,3*numofatoms),float)
# for i in range(len(mass)):
#  mmat[3*i][3*i]=mass[i]**(-0.5)
#  mmat[3*i+1][3*i+1]=mass[i]**(-0.5)
#  mmat[3*i+2][3*i+2]=mass[i]**(-0.5)
# eigvect=dot(eigvect,mmat)
# for i in range(3*numofatoms):
#  eigvect[i]/=sum((eigvect[i]**2)**0.5)


# f=open('vibModes.aclimax','w')
#rrr= "BEGIN ATOMS"+'\r\n'+'\r\n'
# f.write(rrr)
#
# for i in range(numofatoms):
#  amass=mass[i]/amutokg
#  atag=attags[i]
#  asigma=sigmas[attags[i]]
#  apte=pte[attags[i]]
#  # take care of D
#  # Lubo recommended the value of 2.05 but the aCLIMAX default is 7.6!!!
#  if (atag=='H' and amass>1.9):
#    atag='D'
#    asigma=2.05
#  rrr= str(i+1)+" "+atag+" "+str(apte)+" "+'%.3f'%(round(amass,6))
#  rrr+=" "+'%.6f'%(round(coords_c[i][0],6))
#  rrr+=" "+'%.6f'%(round(coords_c[i][1],6))
#  rrr+=" "+'%.6f'%(round(coords_c[i][2],6))
#  rrr+=" "+'%.3f'%(round(asigma,6))
#  rrr+='\r\n'
#  f.write(rrr)
#rrr= '\r\n'+"END ATOMS"+'\r\n'
# f.write(rrr)
#
#
#rrr= '\r\n'+"BEGIN FREQUENCIES"+'\r\n'
# f.write(rrr)
#
# count=0
# for i in range(len(eigval)):
#  # write only real frequencies in cm-1!!!
#  if eigval[i]>0.:
#    ###if (eigval[i]**0.5)/cl/(2*math.pi)*1e-02>1e-1:
#      count+=1
#      rrr='\r\n'+str(count)+" "
#      rrr+='%.2f'%round((eigval[i]**0.5)/cl/(2*math.pi)*1e-02,2)+'\r\n'+'\r\n'
#      f.write(rrr)
#      for j in range(numofatoms):
#        apte=pte[attags[j]]
#        indx=3*j
#        rrr= str(j+1)+" "+attags[j]+" "+str(apte)+" "
#        rrr+= '%.6f'%(round(eigvect[i][indx],6))
#        rrr+=" "+'%.6f'%(round(eigvect[i][indx+1],6))
#        rrr+=" "+'%.6f'%(round(eigvect[i][indx+2],6))+'\r\n'
#        f.write(rrr)
#
#
#rrr= '\r\n'+"END FREQUENCIES"+'\r\n'
# f.write(rrr)
#
#rrr= '\r\n'+"BEGIN PARAMETERS"+'\r\n'
# f.write(rrr)
#rrr= '\r\n'+"END PARAMETERS"+'\r\n'
# f.write(rrr)
#
# f.close()
#
# c from now on the eigvect contains
# c again eigenvectors of dynamical matrix
# eigvect=1*eigvect_
