# phet and i, R is distance variable for Rvals memory 1000 mb molecule dimer { -1 1 C1 -1.845095 0.000000 0.000000 H2 -2.940613 0.000000 0.000000 C3 -1.113628 1.209982 0.000000 H4 -1.651511 2.168601 0.000000 C5 -1.113628 -1.209983 0.000000 H6 -1.651510 -2.168602 0.000000 C7 0.291542 1.215097 0.000000 H8 0.838172 2.167255 0.000000 C9 0.291542 -1.215097 0.000000 H10 0.838172 -2.167254 0.000000 C11 1.086579 0.000000 0.000000 O12 2.372926 0.000000 0.000000 -- -1 1 i O12 R C11 120.0 C9 180.0 } # set the scan variable of R above, ex: distance, manually/automatically assign # ---------------------------------------------------------------- def frange(start, stop, step): ## do not change the section of def frange i = start while i < stop: yield i i += step # ---------------------------------------------------------------- # Rvals=[2.5, 3.0, 4.0] ## manually assign the variable, ex: distance at 2.5, 3.0 and 4.0 anstrom Rvals=[] for j in frange(2.0, 4.0, 0.1): ## automatically assign the variable; the endpoint will not be performed Rvals.append(round(j,2)) ## decimal is 2 # basis set # ---------------------------------------------------------------- ##set basis basis { assign H def2-tzvppd assign Li def2-tzvppd assign Be def2-tzvppd assign C def2-tzvppd assign N def2-tzvppd assign O def2-tzvppd assign F def2-tzvppd assign Na def2-tzvppd assign Mg def2-tzvppd assign P def2-tzvppd assign S def2-tzvppd assign Cl def2-tzvppd assign K def2-tzvppd assign Ca def2-tzvppd assign Br def2-tzvppd assign I def2-tzvppd } set guess sad set scf_type df set maxiter 1000 set mom_start 31 set soscf true set soscf_max_iter 15 set soscf_conv 1.e-4 set freeze_core false ## false for ions ####cp("df-mp2") ## moved to below # get the energy at each variable # ---------------------------------------------------------------------------- # Initialize a blank dictionary of counterpoise corrected energies # (Need this for the syntax below to work) ecp = {} for R in Rvals: dimer.R = R ecp[R] = energy("mp2", bsse_type = "cp") ## for new version #ecp[R] = cp("df-mp2") ## for beta version e= ecp[R] * psi_hartree2kcalmol psi4.print_out("R, E_int [kcal/mol]: %3.1f %10.6f\n" % (R, e)) # print out all the energies # ---------------------------------------------------------------------------- psi4.print_out("\n") psi4.print_out("CP-corrected interaction energies\n\n") psi4.print_out(" R [Ang] E_int [kcal/mol] \n") psi4.print_out("-----------------------------------------------------\n") for R in Rvals: e = ecp[R] * psi_hartree2kcalmol psi4.print_out(" %3.1f %10.6f\n" % (R, e))