diff --git a/bin/SHARC_BAGEL.py b/bin/SHARC_BAGEL.py index bcbf11f..d50c3a2 100755 --- a/bin/SHARC_BAGEL.py +++ b/bin/SHARC_BAGEL.py @@ -1599,20 +1599,20 @@ def readQMin(QMinfilename): sys.exit(44) # PyQuante - QMin['pyquante'] = get_sh2BAGEL_environ(sh2BAGEL, 'pyquante', False, False) - if QMin['pyquante'] is None or not os.path.isdir(QMin['pyquante']): - print('Give path to the PyQuante installation directory in BAGEL.resources!') - sys.exit(45) - if 'PYTHONPATH' in os.environ: - # os.environ['PYTHONPATH']=os.path.join(QMin['pyquante']) + os.pathsep + os.environ['PYTHONPATH'] - - # os.pathsep + QMin['pyquante'] - # print("already pypath") - # print(os.environ['PYTHONPATH']) - # else: - # os.environ['PYTHONPATH']=QMin['pyquante'] - # print("no pypath") - sys.path.append(QMin['pyquante']) + # QMin['pyquante'] = get_sh2BAGEL_environ(sh2BAGEL, 'pyquante', False, False) + # if QMin['pyquante'] is None or not os.path.isdir(QMin['pyquante']): + # print('Give path to the PyQuante installation directory in BAGEL.resources!') + # sys.exit(45) + # if 'PYTHONPATH' in os.environ: + # # os.environ['PYTHONPATH']=os.path.join(QMin['pyquante']) + os.pathsep + os.environ['PYTHONPATH'] + + # # os.pathsep + QMin['pyquante'] + # # print("already pypath") + # # print(os.environ['PYTHONPATH']) + # # else: + # # os.environ['PYTHONPATH']=QMin['pyquante'] + # # print("no pypath") + # sys.path.append(QMin['pyquante']) # memory cannot set memory! #TODO QMin['memory'] = 100 @@ -2334,6 +2334,7 @@ def writeBAGELinput(QMin): # molecule string += ' { \n "title" : "molecule",\n' string += ' "%s" : "%s",\n' % ('dkh', QMin['template']['dkh']) + string += ' "cartesian": "true",\n' # basis if not os.path.isabs(QMin['template']['basis']): @@ -2369,6 +2370,32 @@ def writeBAGELinput(QMin): casscfkeys = ['nact', 'nclosed', 'maxiter'] caspt2keys = ['ms', 'xms', 'sssr', 'shift', 'shift_imag', 'orthogonal_basis', 'maxiter'] + # we must do a force-less calculation first, otherwise we get garbage Molden files + string += ' {\n "title" : "casscf",\n' + for key in casscfkeys: + string += ' "%s" : %s,\n' % (key, QMin['template'][key]) + string += ' "nspin" : %i,\n' % (gsmult - 1) + string += ' "print_thresh" : 1e-10,\n' + string += ' "charge" : %i,\n' % (charge) + string += ' "nstate" : %i,\n' % (QMin['template']['nstate'][gsmult - 1]) + string = string[:-2] + '\n },\n\n' + if QMin['template']['method'] == 'caspt2': + string += ' { \n ' + string += ' "title" : "smith",\n' + string += ' "method" : "caspt2",\n' + for key in caspt2keys: + string += ' "%s" : "%s",\n' % (key, QMin['template'][key]) + string = string[:-2] + '\n },\n\n' + + # save_ref + string += ' { \n "title" : "save_ref",\n' + string += ' "file" : "%s"\n },\n\n' % (os.path.join('./archive')) + + # write molden file + string += ' { \n "title" : "print",\n' + string += ' "file" : "%s",\n' % (os.path.join('./orbitals.molden.%i' % QMin['IJOB'])) + string += ' "orbitals" : "true"\n },\n\n' + # forces if dograd or donac: string += ' { \n "title" : "forces",\n' @@ -2397,32 +2424,24 @@ def writeBAGELinput(QMin): string += ' "nstate" : %i,\n' % (QMin['template']['nstate'][gsmult - 1]) string += ' "dipole" : "%s",\n' % (QMin['template']['dipole']) string = string[:-2] + '\n } ]\n },\n\n' - else: - string += ' {\n "title" : "casscf",\n' - for key in casscfkeys: - string += ' "%s" : %s,\n' % (key, QMin['template'][key]) - string += ' "nspin" : %i,\n' % (gsmult - 1) - string += ' "print_thresh" : 1e-10,\n' - string += ' "charge" : %i,\n' % (charge) - string += ' "nstate" : %i,\n' % (QMin['template']['nstate'][gsmult - 1]) - string = string[:-2] + '\n },\n\n' - if QMin['template']['method'] == 'caspt2': - string += ' { \n ' - string += ' "title" : "smith",\n' - string += ' "method" : "caspt2",\n' - for key in caspt2keys: - string += ' "%s" : "%s",\n' % (key, QMin['template'][key]) - string = string[:-2] + '\n },\n\n' + #else: + # string += ' {\n "title" : "casscf",\n' + # for key in casscfkeys: + # string += ' "%s" : %s,\n' % (key, QMin['template'][key]) + # string += ' "nspin" : %i,\n' % (gsmult - 1) + # string += ' "print_thresh" : 1e-10,\n' + # string += ' "charge" : %i,\n' % (charge) + # string += ' "nstate" : %i,\n' % (QMin['template']['nstate'][gsmult - 1]) + # string = string[:-2] + '\n },\n\n' + # if QMin['template']['method'] == 'caspt2': + # string += ' { \n ' + # string += ' "title" : "smith",\n' + # string += ' "method" : "caspt2",\n' + # for key in caspt2keys: + # string += ' "%s" : "%s",\n' % (key, QMin['template'][key]) + # string = string[:-2] + '\n },\n\n' - # save_ref - string += ' { \n "title" : "save_ref",\n' - string += ' "file" : "%s"\n },\n\n' % (os.path.join('./archive')) - - # write molden file - string += ' { \n "title" : "print",\n' - string += ' "file" : "%s",\n' % (os.path.join('./orbitals.molden.%i' % QMin['IJOB'])) - string += ' "orbitals" : "true"\n },\n\n' string = string[:-3] + '\n\n' string += ']}' @@ -2697,6 +2716,8 @@ def make_mos_from_Molden(moldenfile, QMin): NAOs[1] += shells[s[0]][mode[s[0]]] aos.append(s[0]) + #print(NAOs) + #print(aos) for iline, line in enumerate(data): if '[mo]' in line.lower(): @@ -2716,59 +2737,90 @@ def make_mos_from_Molden(moldenfile, QMin): MO_A[imo][iao] = float(line.split()[1]) if any([1 == mode[i] for i in aos]): - print(' .. transforming to cartesian AO basis') - - - d34 = math.sqrt(3. / 4.) - f65 = math.sqrt(6. / 5.) - f38 = math.sqrt(3. / 8.) - f58 = math.sqrt(5. / 8.) - f98 = math.sqrt(9. / 8.) - f920 = math.sqrt(9. / 20.) - f340 = math.sqrt(3. / 40.) - transform = { - 'd': [[(0, -0.5), (3, d34)], - [(0, -0.5), (3, -d34)], - [(0, 1.0)], - [(4, 1.0)], - [(1, 1.0)], - [(2, 1.0)]], - 'f': [[(1, -f38), (5, f58)], - [(2, -f38), (6, -f58)], - [(0, 1.0)], - [(1, -f340), (5, -f98)], - [(2, -f340), (6, f98)], - [(0, -f920), (3, d34)], - [(1, f65)], - [(2, f65)], - [(0, -f920), (3, -d34)], - [(4, 1.0)]], - 'g': None # TODO: add and check for bagel + print('no support for spherical basis sets in overlaps') + sys.exit(1) + + # reorder AOs from Molden order to pyscf order + # https://www.theochem.ru.nl/molden/molden_format.html + # https://pyscf.org/user/gto.html#ordering-of-basis-function + reorders = { + 's': {0:0}, + 'p': {0:0, 1:1, 2:2}, + 'd': {0:0, 1:3, 2:5, 3:1, 4:2, 5:4}, + 'f': {0:0, 1:6, 2:9, 3:3, 4:1, 5:2, 6:4, 7:8, 8:7, 9:4} } - - tMO_A = [] + factors = { + #'d': {1: math.sqrt(3.), 2: math.sqrt(3.), 4: math.sqrt(3.)}, + } + MO_A_reorder = [] for imo in range(NMO_A): - newmo = [] - iao = 0 - for shell in aos: - n_cart = shells[shell][0] - if mode[shell] == 0: - for i in range(n_cart): - newmo.append(MO_A[imo][iao]) - iao += 1 - else: - n_sph = shells[shell][1] - new = [0. for i in range(n_cart)] - for i in range(n_cart): - for j, c in transform[shell][i]: - new[i] += c * MO_A[imo][iao + j] - newmo.extend(new) - iao += n_sph - tMO_A.append(newmo) - - - MO_A = tMO_A - NAO = NAOs[0] + mo = [0. for i in range(NAO)] + ishift=0 + for ishell,s in enumerate(aos): + for i in range(shells[s][mode[s[0]]]): + iao_old = ishift + i + iao_new = ishift + reorders[s][i] + #print(ishell,s,i,iao_old,'->',iao_new)i + if s in factors and i in factors[s]: + fact=factors[s][i] + else: + fact=1. + mo[iao_new] = MO_A[imo][iao_old]*math.sqrt(fact) + ishift += shells[s][mode[s[0]]] + MO_A_reorder.append(mo) + MO_A = MO_A_reorder + + + #d34 = math.sqrt(3. / 4.) + #f65 = math.sqrt(6. / 5.) + #f38 = math.sqrt(3. / 8.) + #f58 = math.sqrt(5. / 8.) + #f98 = math.sqrt(9. / 8.) + #f920 = math.sqrt(9. / 20.) + #f340 = math.sqrt(3. / 40.) + #transform = { + # 'd': [[(0, -0.5), (3, d34)], + # [(0, -0.5), (3, -d34)], + # [(0, 1.0)], + # [(4, 1.0)], + # [(1, 1.0)], + # [(2, 1.0)]], + # 'f': [[(1, -f38), (5, f58)], + # [(2, -f38), (6, -f58)], + # [(0, 1.0)], + # [(1, -f340), (5, -f98)], + # [(2, -f340), (6, f98)], + # [(0, -f920), (3, d34)], + # [(1, f65)], + # [(2, f65)], + # [(0, -f920), (3, -d34)], + # [(4, 1.0)]], + # 'g': None # TODO: add and check for bagel + #} + + #tMO_A = [] + #for imo in range(NMO_A): + # newmo = [] + # iao = 0 + # for shell in aos: + # n_cart = shells[shell][0] + # if mode[shell] == 0: + # for i in range(n_cart): + # newmo.append(MO_A[imo][iao]) + # iao += 1 + # else: + # n_sph = shells[shell][1] + # new = [0. for i in range(n_cart)] + # for i in range(n_cart): + # for j, c in transform[shell][i]: + # new[i] += c * MO_A[imo][iao + j] + # newmo.extend(new) + # iao += n_sph + # tMO_A.append(newmo) + + + #MO_A = tMO_A + #NAO = NAOs[0] # handle frozen core @@ -3259,9 +3311,15 @@ def get_smat_from_Molden(file1, file2=''): mol2, mo_energy, mo_coeff, mo_occ, irrep_labels, spins = load(file2) mol = mole.conc_mol(mol, mol2) - S = mol.intor('int1e_ovlp').tolist() + S = mol.intor('int1e_ovlp') #.tolist() + + # normalize with sqrt(diag) + for i in range(len(S)): + n=math.sqrt(S[i,i]) + S[i,:] /= n + S[:,i] /= n - return len(S), S + return len(S), S.tolist() # ======================================================================= # @@ -3603,7 +3661,7 @@ def getenergy(corefile, ijob, QMin): i = int(line.split()[5]) energies[(imult, i + (gsmult == imult))] = float(line.split()[6]) - else: + elif not mscaspt2 and not caspt2: if '* ci vector' in line and ' 0,' in line: for i in range(nstates): if i == 0: @@ -3934,7 +3992,7 @@ def main(): # Remove Scratchfiles from SCRATCHDIR if not DEBUG: - cleandir(QMin['scratchdir']) + #cleandir(QMin['scratchdir']) if 'cleanup' in QMin: cleandir(QMin['savedir']) diff --git a/bin/SHARC_MOLCAS.py b/bin/SHARC_MOLCAS.py index 92a552f..e0550ba 100755 --- a/bin/SHARC_MOLCAS.py +++ b/bin/SHARC_MOLCAS.py @@ -1042,7 +1042,8 @@ def getcidm(out, mult, state1, state2, pol, version): stop2string = 'Special properties section' statesstring = 'Nr of states:' #matrixstring=' PROPERTY: MLTPL 1 COMPONENT: %i' % (pol+1) - matrixstring = 'PROPERTY: MLTPL 1 COMPONENT:* %i' % (pol + 1) + #matrixstring = 'PROPERTY: MLTPL 1 COMPONENT:* %i' % (pol + 1) + matrixstring = 'PROPERTY: MLTPL 1 COMPONENT:' # first, find the correct RASSI output section for the given multiplicity module = False @@ -1077,7 +1078,7 @@ def getcidm(out, mult, state1, state2, pol, version): stateshift = nstates // 2 else: stateshift = 0 - if matrixstring in line: + if matrixstring in line and int(line.split()[-1]) == pol+1: block = (stateshift + state2 - 1) // 4 rowshift = 3 + stateshift + state1 + (6 + nstates) * block colshift = 1 + (stateshift + state2 - 1) % 4 @@ -2240,12 +2241,14 @@ def readQMin(QMinfilename): driver = get_sh2cas_environ(sh2cas, 'driver', crucial=False) if driver == '': - driver = os.path.join(QMin['molcas'], 'bin', 'pymolcas') + driver = os.path.join(QMin['molcas'], 'pymolcas') if not os.path.isfile(driver): - driver = os.path.join(QMin['molcas'], 'bin', 'molcas.exe') + driver = os.path.join(QMin['molcas'],'bin', 'pymolcas') if not os.path.isfile(driver): - print('No driver (pymolcas or molcas.exe) found in $MOLCAS/bin. Please add the path to the driver via the "driver" keyword.') - sys.exit(52) + driver = os.path.join(QMin['molcas'], 'bin', 'molcas.exe') + if not os.path.isfile(driver): + print('No driver (pymolcas or molcas.exe) found in $MOLCAS or $MOLCAS/bin. Please add the path to the driver via the "driver" keyword.') + sys.exit(52) QMin['driver'] = driver QMin['tinker'] = get_sh2cas_environ(sh2cas, 'tinker', crucial=False) diff --git a/bin/create_LVCparam.py b/bin/create_LVCparam.py index 88c61b2..1f42f4a 100755 --- a/bin/create_LVCparam.py +++ b/bin/create_LVCparam.py @@ -286,22 +286,32 @@ def LVC_complex_mat(header, mat, deldiag=False, oformat=' % .7e'): # ======================================================================= # -def loewdin_orthonormalization(A): +def loewdin_orthonormalization(A, debug=False): ''' returns loewdin orthonormalized matrix ''' # S = A^T * A S = np.dot(A.T, A) + if debug: + print(S) # S^d = U^T * S * U S_diag_only, U = np.linalg.eigh(S) + if debug: + print('test') + print(U) + print(U.shape) + print(U[0].shape) # calculate the inverse sqrt of the diagonal matrix S_diag_only_inverse_sqrt = [1. / (float(d) ** 0.5) for d in S_diag_only] S_diag_inverse_sqrt = np.diag(S_diag_only_inverse_sqrt) # calculate inverse sqrt of S + if debug: + printmatrix(S_diag_inverse_sqrt) + printmatrix(U.T,'bla') S_inverse_sqrt = np.dot(np.dot(U, S_diag_inverse_sqrt), U.T) # calculate loewdin orthonormalized matrix @@ -312,6 +322,8 @@ def loewdin_orthonormalization(A): length = len(A_lo) A_lon = np.zeros((length, length), dtype=complex) + if debug: + printmatrix(A_lo,'A_lo before normalization') for i in range(length): norm_of_col = np.linalg.norm(A_lo[i]) A_lon[i] = [e / (norm_of_col ** 0.5) for e in A_lo[i]][0] @@ -364,7 +376,7 @@ def partition_matrix(matrix, multiplicity, states): # ======================================================================= # -def phase_correction(matrix): +def phase_correction(matrix, debug=False): U = np.array(matrix).real #print(U) det = np.linalg.det(U) @@ -439,7 +451,7 @@ def check_overlap_diagonal(matrix, states, normal_mode, displacement, ignore_pro # ======================================================================= # -def calculate_W_dQi(H, S, e_ref, normal_mode, displ): +def calculate_W_dQi(H, S, e_ref, normal_mode, displ, debug=False): ''' Calculates the displacement matrix ''' @@ -449,11 +461,22 @@ def calculate_W_dQi(H, S, e_ref, normal_mode, displ): # do phase correction if necessary # if any([x for x in np.diag(S) if x < 0]): + if debug: + print(normal_mode, displ) + printmatrix(S,'before phase 1') S = phase_correction(S) + if debug: + printmatrix(S,'after phase 1') # do loewdin orthonorm. on overlap matrix - U = loewdin_orthonormalization(np.matrix(S)) + if debug: + printmatrix(S,'before Loewdin') + U = loewdin_orthonormalization(np.matrix(S), debug=debug) + if debug: + printmatrix(U,'after Loewdin') U = np.array(phase_correction(U)) + if debug: + printmatrix(U,'after phase 2') return np.dot(np.dot(U, H), U.T) # TODO: or transposed the other way around? # return np.dot(np.dot(U.T, H), U) @@ -463,7 +486,9 @@ def calculate_W_dQi(H, S, e_ref, normal_mode, displ): def printmatrix(M,title=''): s='%s\n' % title for i in M: + print(i) for j in i: + print(j) s+='%6.3f ' % j.real s+='\n' print(s) @@ -617,7 +642,7 @@ def write_LVC_template(INFOS): INFOS['problematic_mults'] = check_overlap_diagonal(pos_S, INFOS['states'], normal_mode, 'p', INFOS['ignore_problematic_states']) # calculate displacement matrix - pos_W_dQi = calculate_W_dQi(pos_H, pos_S, e_ref, normal_mode, 'p') + pos_W_dQi = calculate_W_dQi(pos_H, pos_S, e_ref, normal_mode, 'p', debug=INFOS['debug']) print('Mode %s positive' % normal_mode) # printmatrix(pos_H, 'Hpos') # printmatrix(pos_S, 'Spos') @@ -641,7 +666,7 @@ def write_LVC_template(INFOS): INFOS['problematic_mults'].update(check_overlap_diagonal(neg_S, INFOS['states'], normal_mode, 'n', INFOS['ignore_problematic_states'])) # calculate displacement matrix - neg_W_dQi = calculate_W_dQi(neg_H, neg_S, e_ref, normal_mode, 'n') + neg_W_dQi = calculate_W_dQi(neg_H, neg_S, e_ref, normal_mode, 'n', debug=INFOS['debug']) print('Mode %s negative' % normal_mode) # printmatrix(neg_H, 'Hneg') # printmatrix(neg_S, 'Sneg') @@ -735,6 +760,10 @@ def main(): usage = '''python %s''' % (script_name) parser = OptionParser(usage=usage, description='') + parser.add_option("-d", "--debug", + action="store_true", dest="debug", default=False, + help="Print all involved matrices and quantities for debugging") + (options, args) = parser.parse_args() displaywelcome() @@ -750,6 +779,7 @@ def main(): # set manually for old calcs # INFOS['ignore_problematic_states'] = True + INFOS['debug'] = options.debug # write LVC.template write_LVC_template(INFOS) diff --git a/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF/run.sh b/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF/run.sh new file mode 100755 index 0000000..58a0c91 --- /dev/null +++ b/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF/run.sh @@ -0,0 +1,20 @@ +#/bin/bash + +# $PRIMARYDIR and SCRADIR are set +COPY_DIR=$SCRADIR/TRAJ +PRIMARYDIR=$(pwd) + +mkdir -p $COPY_DIR +cp -r $PRIMARYDIR/* $COPY_DIR +cd $COPY_DIR + +$SHARC/sharc.x input +err=$? + +cp $COPY_DIR/output.* $PRIMARYDIR +if [ ! $err == 0 ]; +then + cp $COPY_DIR/QM/* $PRIMARYDIR/QM/ +fi +rm -r $COPY_DIR +exit $err diff --git a/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF_with_TDH_scheme/run.sh b/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF_with_TDH_scheme/run.sh new file mode 100755 index 0000000..58a0c91 --- /dev/null +++ b/tests/INPUT/MOLCAS_tCSDM_SA-CASSCF_with_TDH_scheme/run.sh @@ -0,0 +1,20 @@ +#/bin/bash + +# $PRIMARYDIR and SCRADIR are set +COPY_DIR=$SCRADIR/TRAJ +PRIMARYDIR=$(pwd) + +mkdir -p $COPY_DIR +cp -r $PRIMARYDIR/* $COPY_DIR +cd $COPY_DIR + +$SHARC/sharc.x input +err=$? + +cp $COPY_DIR/output.* $PRIMARYDIR +if [ ! $err == 0 ]; +then + cp $COPY_DIR/QM/* $PRIMARYDIR/QM/ +fi +rm -r $COPY_DIR +exit $err