hdiff output

r31586/path_test_parser.py 2016-12-05 19:30:11.234316326 +0000 r31585/path_test_parser.py 2016-12-05 19:30:13.118340733 +0000
  1: import sys  1: import sys
  2: import numpy as np 
  3:  
  4: def get_Is(config): 
  5:     Q = np.array(config).reshape(-1,3) 
  6:     CoM = np.average(Q,axis=0) 
  7:     Q-=CoM 
  8:      
  9:     I = np.zeros((3,3))  # non-mass-weighted Inertia tensor 
 10:     for atom in Q: 
 11:         I[0][0]+=atom[1]*atom[1]+atom[2]*atom[2] 
 12:         I[1][1]+=atom[0]*atom[0]+atom[2]*atom[2] 
 13:         I[2][2]+=atom[1]*atom[1]+atom[0]*atom[0] 
 14:         I[0][1]-=atom[0]*atom[1] 
 15:         I[0][2]-=atom[0]*atom[2] 
 16:         I[1][2]-=atom[1]*atom[2] 
 17:     I[1][0] = I[0][1] 
 18:     I[0][2] = I[2][0] 
 19:     I[1][2] = I[2][1] 
 20:  
 21:     return np.linalg.eigvalsh(I)  # The eigenvalues of I 
 22:  
 23:   2: 
 24: def read_pathinfo(fname, nconfigs):  3: def read_pathinfo(fname, nconfigs):
 25:     fin = open(fname,'r')  4:     fin = open(fname,'r')
 26:   5: 
 27:     energies = []  6:     energies = []
 28:     symmetries = []  7:     symmetries = []
 29:     inertias = []  8:     configs = []
 30:     thisconfig = []  9:     thisconfig = []
 31:  10: 
 32:     sline = fin.readline().split() # Read the first energy line 11:     sline = fin.readline().split() # Read the first energy line
 33:     for i in xrange(nconfigs): 12:     for i in xrange(nconfigs):
 34:         if len(sline)!=1: 13:         if len(sline)!=1:
 35:             raise ValueError("Energy line appears to have multiple entries. Check nconfigs = "+str(nconfigs)) 14:             raise ValueError("Energy line appears to have multiple entries. Check nconfigs = "+nconfigs)
 36:         energies.append(float(sline[0]))   # Save the energy of this configuration 15:         energies.append(float(sline[0]))   # Save the energy of this configuration
 37:         symmetries.append(fin.readline())  # Read point group information for this configuration 16:         symmetries.append(fin.readline())  # Read point group information for this configuration
 38:         sline = fin.readline().split()     # Read the first coordinate/frequency line 17:         sline = fin.readline().split()     # Read the first coordinate/frequency line
 39:         while len(sline)==3: 18:         while len(sline)==3:
 40:             thisconfig.append(map(float, sline)) 19:             thisconfig.append(map(float, sline))
 41:             sline = fin.readline().split() # Read the next coordinate/frequency line. 20:             sline = fin.readline().split() # Read the next coordinate/frequency line.
 42:                                            # When this readline() encounters the start of the next configuration, 21:                                            # When this readline() encounters the start of the next configuration,
 43:                                            # we exit the loop with sline containing the next energy line, ready for 22:                                            # we exit the loop with sline containing the next energy line, ready for
 44:                                            # the next cycle of the loop over i. 23:                                            # the next cycle of the loop over i.
 45:         Is = get_Is(thisconfig) 24:         configs.append(thisconfig)  # Save the current configuration
 46:         inertias.append(Is)                # Save principle moments of inertia for the current configuration 25:         thisconfig=[]               # and empty ready to receive the next
 47:         thisconfig=[]                      # and empty ready to receive the next 
 48:  26: 
 49:     fin.close() 27:     fin.close()
 50:  28: 
 51:     return energies, symmetries, inertias 29:     return energies, symmetries, configs
 52:  
 53: def check_ordering(E1, E2, tol): 
 54:     nconfigs = len(E1) 
 55:     mapping = range(nconfigs) 
 56:     for i in xrange(nconfigs): 
 57:         if mapping[i]==i and abs((E1[i]-E2[i])/E1[i]) > tol: 
 58:             for j in xrange(i+1, nconfigs): 
 59:                 if abs((E1[i]-E2[j])/E1[i]) < tol: 
 60:                     mapping[i] = j 
 61:                     mapping[j] = i 
 62:                     break 
 63:             if mapping[i]==i: 
 64:                 print "Failed to find a stationary point with the same energy as ", i 
 65:                 print "Mapping points together is not possible." 
 66:                 return None  # Failure: there is no configuration which matches configuration i. 
 67:     return mapping 
 68:  
 69:  30: 
 70: tol = 1.0e-4 31: tol = 1.0e-4
 71: fname1 = 'path.info' 32: fname1 = 'path.info'
 72: fname2 = 'expected_output/path.info' 33: fname2 = 'expected_output/path.info'
 73: print "Comparing ", fname1, "with ", fname2 34: print "Comparing ", fname1, "with ", fname2
 74: success = True 35: success = True
 75:  36: 
 76: # We are assuming only one triple per path.info file. This is a pretty safe bet, since the PATH run should only ever produce one triple. 37: # We are assuming only one triple per path.info file. This is a pretty safe bet, since the PATH run should only ever produce one triple.
 77: E1, sym1, I1 = read_pathinfo(fname1, 3) 38: E1, sym1, configs1 = read_pathinfo(fname1, 3)
 78: E2, sym2, I2 = read_pathinfo(fname2, 3) 39: E2, sym2, configs2 = read_pathinfo(fname2, 3)
 79:      40:     
  41: natoms = len(configs1[0]) # Since we have no frequencies printed here, the number of coordinate lines is the same as the number
  42:                           # of atoms.
  43: print "natoms ", natoms
  44: if len(configs2) != len(configs1):
  45:     print "The path.info files appear to belong to different systems."
  46:     success = False
  47: 
 80: nconfigs = len(E1) 48: nconfigs = len(E1)
 81: print "nconfigs", nconfigs 49: print "nconfigs", nconfigs # This should always be 3 for PATH tests
 82: if len(E2)!=len(E1): 50: if len(E2)!=len(E1):
 83:     print "The path.info have different numbers of entries." 51:     print "The path.info have different numbers of entries."
 84:     success = False 52:     success = False
 85:  53: 
 86: if success: 54: if success:
 87:     Efailures = 0 55:     Efailures = 0
 88:     symfailures = 0 56:     symfailures = 0
 89:     Ifailures = 0 57:     coordsfailures = 0
 90:  58:     coordsfailconfigs = {}
 91:     mapping = check_ordering(E1, E2, tol) 
 92:  
 93:     if mapping is not None and mapping!=range(nconfigs):  # if it is None, then we will report the failure in the normal way 
 94:         print "Mapping stationary points between the files according to: ", mapping 
 95:         E2 = [ E2[i] for i in mapping] 
 96:         sym2 = [ sym2[i] for i in mapping] 
 97:         I2 = [ I2[i] for i in mapping]        
 98:  59: 
 99:     for i in xrange(nconfigs): 60:     for i in xrange(nconfigs):
100:         if abs((E1[i]-E2[i])/E1[i]) > tol: 61:         if abs((E1[i]-E2[i])/E1[i]) > tol:
101:             Efailures +=1 62:             Efailures +=1
102:             print "Different energies for structure "+str(i)+':', E1[i], E2[i] 63:             print "Different energies for structure "+str(i)+':', E1[i], E2[i]
103:             success = False 64:             success = False
104:  65: 
105:     for i in xrange(nconfigs): 66:     for i in xrange(nconfigs):
106:         if sym1[i]!=sym2[i]: 67:         if sym1[i]!=sym2[i]:
107:             symfailures+=1 68:             symfailures+=1
108:             print "Symmetries differ for structure "+str(i)+':', sym1[i][:-2], sym2[i][:-2] 69:             print "Symmetries differ for structure "+str(i)+':', sym1[i][:-2], sym2[i][:-2]
109:             success = False 70:             success = False
110:  71: 
111:     for i in xrange(nconfigs): 72:     for i in xrange(nconfigs):
112:         if (abs((I1[i][0]-I2[i][0])/I1[i][0])>tol) or (abs((I1[i][1]-I2[i][1])/I1[i][1])>tol) or (abs((I1[i][2]-I2[i][2])/I1[i][2])>tol): 73:         coordsfailconfigs[i] = 0
113:             Ifailures+=1 74:         for j in xrange(natoms):
114:             print "Different moments of inertia for structure "+str(i)+':' 75:             for k in xrange(3):
115:             print I1[i] 76:                 if abs((configs1[i][j][k]-configs2[i][j][k])/configs1[i][j][k]) > tol:
116:             print I2[i] 77:                     coordsfailures+=1
  78:                     coordsfailconfigs[i]+=1
  79: 
  80:         if coordsfailconfigs[i]!=0:
  81:             print coordsfailconfigs[i], "coordinates changed for structure ", i
117:             success = False 82:             success = False
118:  83: 
119: else:  # If the first safety checks failed 84: else:  # If the first safety checks failed
120:     print "path.info files differed dramatically: comparison not possible." 85:     print "path.info files differed dramatically: comparison not possible."
121:     sys.exit(1) 86:     sys.exit(1)
122:  87: 
123:  88: 
124: if success: 89: if success:
125:     print "All tests completed successfully." 90:     print "All tests completed successfully."
126:     sys.exit(0) 91:     sys.exit(0)
127: else: 92: else:
128:     print "Number of stationary point energies which differed between the files:", Efailures 93:     print "Number of stationary point energies which differed between the files:", Efailures
129:     print "Number of stationary point symmetries which differed between the files:", symfailures 94:     print "Number of stationary point symmetries which differed between the files:", symfailures
130:     print "Number of stationary points where the moments of inertia differed between the files:", Ifailures 95:     print "Number of coordinates which differed between the files:", coordsfailures
131:     sys.exit(1) 96:     sys.exit(1)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0