FIXME: virial theorem checking, but must calculate averages. So it's …

…to be fixed.
cosurgi · Jul 1, 2015 · 2b63787 · 2b63787
1 parent f0ec392
commit 2b63787
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 0 deletions.
diff --git a/pkg/qm/SchrodingerPropagator.cpp b/pkg/qm/SchrodingerPropagator.cpp
@@ -103,6 +103,12 @@ o̲n̲ ̲e̲a̲c̲h̲ ̲c̲a̲l̲l̲!̲ ̲I̲ ̲n̲e̲e̲d̲ ̲s̲o̲m̲e̲ ̲d
 
 	//if(allPotentials.size() == 0)	return 0;
 	//NDimTable<Complexr>& Vpsi=(*(allPotentials.begin()))->psiGlobalTable;
+
+//	for(int i=0;i<Vpsi.dim()[0];i++)
+//		std::cerr << i << " " << Vpsi.at(i,Vpsi.dim()[0]/2) << "\n";
+//	for(int i=0;i<Vpsi.dim()[0];i++)
+//		std::cerr << i << " " << Vpsi.at(i) << "\n";
+
 	return std::move(Vpsi);
 };
 
@@ -133,6 +139,20 @@ boost::shared_ptr<QMStateDiscreteGlobal> SchrodingerKosloffPropagator::get_full_
 	return *allPsiGlobals.begin();
 };
 
+// Virial theorem: 2*<E_kin> = -<E_pot>
+
+void SchrodingerKosloffPropagator::virialTheorem_Grid_check()
+{
+	if(virialCheck) {
+		Real E_pot = abs(eMin()); // FIXME - tu powinna być liczona średnia! A nie MAX.
+		Real E_kin = abs(eMax()); // FIXME - tu powinna być liczona średnia! A nie MAX.
+		if(2*E_kin > E_pot) {
+			std::cerr << "ERROR: too big potential energy for this gridsize\n";
+			//exit(0);
+		};
+	};
+};
+
 Real SchrodingerKosloffPropagator::eMin()
 {
 	NDimTable<Complexr> VGlobal( get_full_potentialInteractionGlobal_psiGlobalTable() );
@@ -214,6 +234,7 @@ void SchrodingerKosloffPropagator::calc_Hnorm_psi(const NDimTable<Complexr>& psi
 
 void SchrodingerKosloffPropagator::action()
 {
+	virialTheorem_Grid_check(); // FIXME - to powinno być chyba zależne od potencjału...
 	timeLimit.readWallClock();
 	Real R   = calcKosloffR(scene->dt); // FIXME -  that's duplicate here, depends on dt !!
 	Real G   = calcKosloffG(scene->dt); // FIXME -  that's duplicate here, depends on dt !!

diff --git a/pkg/qm/SchrodingerPropagator.hpp b/pkg/qm/SchrodingerPropagator.hpp
@@ -93,6 +93,7 @@ found in [TalEzer1984]_"
 			, // attributes, public variables
 			((Real    ,hbar,1               ,,"Planck's constant $h$ divided by $2\\pi$"))
 			((int     ,steps     ,-1     ,,"Override automatic selection of number of steps in Chebyshev expansion."))
+			((bool    ,virialCheck,false ,,"Check energies using virial theorem (Coulomb potential ONLY - FIXME!!!!!!!!)."))
 			, // constructor
 			, // python bindings
 			.def("eMin"  ,&SchrodingerKosloffPropagator::eMin  ,"Get minimum energy.")
@@ -110,6 +111,7 @@ found in [TalEzer1984]_"
 		// FIXME są różne typy, to jest podejrzane. Może w QMIPhys wystarczy trzymać NDimTable, a nie całe QMStateDiscreteGlobal ?
 		NDimTable<Complexr>                      get_full_potentialInteractionGlobal_psiGlobalTable();
 		boost::shared_ptr<QMStateDiscreteGlobal> get_full_psiGlobal__________________psiGlobalTable();
+		void virialTheorem_Grid_check();
 };
 REGISTER_SERIALIZABLE(SchrodingerKosloffPropagator);