incorporated PR comments

Author: msnik1999

external/upstream/fetch_mrcpp.cmake

@@ -39,7 +39,7 @@ else()
     GIT_REPOSITORY
       https://github.com/MRChemSoft/mrcpp.git
     GIT_TAG
-      720133372c9717134c5a01e963cb9804a1e8c36e
+      f9e50bae7c28ece70fa0c909e4be7170d406543a
   )
 
   FetchContent_GetProperties(mrcpp_sources)

src/driver.cpp

@@ -115,6 +115,7 @@ bool guess_energy(const json &input, Molecule &mol, FockBuilder &F);
 void write_orbitals(const json &input, Molecule &mol);
 void calc_properties(const json &input, Molecule &mol);
 void plot_quantities(const json &input, Molecule &mol);
+bool _useExchange = false;
 } // namespace scf
 
 namespace rsp {
@@ -255,7 +256,11 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
     OrbitalVector Phi_mom;
     if (scf::guess_orbitals(json_guess, json_occ, mol)) {
         if (json_scf.contains("scf_solver")) {
-            if (json_scf["scf_solver"]["deltascf_method"] == "IMOM" || json_scf["scf_solver"]["deltascf_method"] == "MOM") Phi_mom = orbital::deep_copy(mol.getOrbitals());
+            if (json_scf["scf_solver"]["deltascf_method"] == "IMOM" || json_scf["scf_solver"]["deltascf_method"] == "MOM") {
+                if (_useExchange)
+                    MSG_ABORT("Running DeltaSCF calculations with exact exchange is currently not supported!");
+                Phi_mom = orbital::deep_copy(mol.getOrbitals());
+            }
         }
         scf::guess_energy(json_guess, mol, F);
         json_out["initial_energy"] = mol.getSCFEnergy().json();
@@ -1275,6 +1280,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
             auto K_p = std::make_shared<ExchangeOperator>(P_p, Phi_p, X_p, Y_p, exchange_prec);
             F.getExchangeOperator() = K_p;
         }
+        scf::_useExchange = true; // determine if exact exchange is used in order to prevent MOM/IMOM calculations
     }
     ///////////////////////////////////////////////////////////
     /////////////////   External Operator   ///////////////////

src/scf_solver/GroundStateSolver.cpp

@@ -324,11 +324,10 @@ json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F, OrbitalVector &P
             if (restricted) {
                 DoubleVector occNew = getNewOccupations(Phi_n, Phi_mom);
                 orbital::set_occupations(Phi_n, occNew);
-                if (plevel > 1) {
-                    orbital::print(Phi_n);
-                    mol.calculateOrbitalPositions();
-                    mol.printOrbitalPositions();
-                }
+                // ToDo: output depending on print level
+                orbital::print(Phi_n);
+                mol.calculateOrbitalPositions();
+                mol.printOrbitalPositions();
             }
             else {
                 // in case of unrestricted calculation, get the new occupation for alpha and beta spins independently
@@ -341,11 +340,10 @@ json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F, OrbitalVector &P
                 DoubleVector occNew(occAlpha.size() + occBeta.size());
                 occNew << occAlpha, occBeta;
                 orbital::set_occupations(Phi_n, occNew);
-                if (plevel > 1) {
-                    orbital::print(Phi_n);
-                    mol.calculateOrbitalPositions();
-                    mol.printOrbitalPositions();
-                }
+                // ToDo: output depending on print level
+                orbital::print(Phi_n);
+                mol.calculateOrbitalPositions();
+                mol.printOrbitalPositions();
             }
         }