2NDforceconstants for VCA

#!/bin/python

# this simple script is to give the combination of two different elements for alloy to calculate the harmonic dynamic properties in almabte
# requirements: sencond order force constants for two different elements and put into current directory
# for example: FORCE_CONSTANTS_2ND_Si and  FORCE_CONSTANTS_2ND_Ge to generate FORCE_CONSTANTS_VCA
# then apply phonopy to get what you want.

import numpy as np
with open('FORCE_CONSTANTS_2ND_Si') as fcfile1:
     with open('FORCE_CONSTANTS_2ND_Ge') as fcfile2:
          idx1 = []
          idx2 = []
          line1 = fcfile1.readline()
          line2 = fcfile2.readline()
          idx1 = [int(x) for x in line1.split()]
          idx2 = [int(y) for y in line2.split()]
          if len(idx1) == 1:
             idx1 = [idx1[0], idx1[0]]
          if len(idx2) == 1:
             idx2 = [idx2[0], idx2[0]]
          force_constants1 = np.zeros((idx1[0], idx1[1], 3, 3), dtype='double')
          force_constants2 = np.zeros((idx2[0], idx2[1], 3, 3), dtype='double')
          force_constants = np.zeros((idx1[0], idx1[1], 3, 3), dtype='double')
          print np.shape(force_constants)
          for i in range(idx1[0]):
              for j in range(idx1[1]):
                  s_i1 = int(fcfile1.readline().split()[0]) - 1
                  s_i2 = int(fcfile2.readline().split()[0]) - 1
                  if s_i1 not in idx1:
                     idx1.append(s_i1)
                  if s_i2 not in idx2:
                     idx2.append(s_i2)
                #print idx1
                  tensor = []
                # here you can change the ratio to what you want
                  for k in range(3):
                      tensor.append([float(x)*0.4 + float(y)*0.6
                                    for x, y in zip(fcfile1.readline().split(), fcfile2.readline().split())])
                  force_constants[i,j] = tensor
             #for l in range(3):
             #    for m in range(3):
             #        print force_constants[i,j,l,m]
fcfile1.close()
fcfile2.close()
indices = np.arange(force_constants.shape[0], dtype='intc')
lines = []
fc_shape = force_constants.shape
print fc_shape
lines.append("%4d %4d" % fc_shape[:2])
for l, s_l in enumerate(indices):
    for m in range(fc_shape[1]):
        lines.append("%d %d" % (s_l + 1, m + 1))
        for vec in force_constants[l][m]:
            lines.append(("%22.15f" * 3) % tuple(vec))
print np.shape(lines), lines[1], type(lines)
with open('FORCE_CONSTANTS_VCA', 'w') as w:
     w.write("\n".join(lines))