Revert singe system interface

- the added code complexity really isn't worth the front end simplifciaton. Going to revert
- back to all lists for mean-field calculations. tutorial mostly done except for a late-binding bug.
- fixed binding issue. Found a few typos including wrong dissipation rate
- documentation and basic cosmetics
- update mean-field tempo example, and check it matches output before changes (difference ~1e-15 in imaginary part can be accounted for just by taking an element of a list e.g. states=[state];state=states[0]!)
- update process example too
- tidy up example with multiple subssytems
- give these examples a realistic hamiltonian
- remove my mean field correlation scripts which shouldn't be in public repo
- rename TempoWithField to MeanFieldTempo
- updated info on PATHs for tests

CONTRIBUTING.md

@@ -204,7 +204,7 @@ and you can pick a specific test to run with:
 $ tox -e py36 say_hi_test.py
 ```
 Here `say_hi_test.py` is the path of the test file *relative to the tests
-directory.*
+directory and this command should be run from the OQuPy base directory!*
 
 #### 7.3 test pylint only
 This checks the code [code style](https://www.python.org/dev/peps/pep-0008/)

examples/mean-field-process.py

@@ -25,18 +25,19 @@
 wc=0.0
 kappa=0.1
 end_time=1
-def field_eom(t, state, field):
-    sx_exp = np.matmul(Sx, state).trace().real
+def field_eom(t, states, field):
+    sx_exp = np.matmul(Sx, states[0]).trace().real
     return -(1j*wc+kappa)*field - 1j*gn*sx_exp
-def H(t, field):
-    #print(t, field)
-    return 2.0 * gn * np.abs(field) * Sx
+def H_MF(t, field):
+    return 2.0 * gn * np.real(field) * Sx
 
-system = oqupy.TimeDependentSystemWithField(H,
-        field_eom,
+system = oqupy.TimeDependentSystemWithField(H_MF,
         #gammas = [gam_down, gam_up],
         #lindblad_operators = [oqupy.operators.sigma("-"), tempo.operators.sigma("+")]
         )
+
+mean_field_system = oqupy.MeanFieldSystem([system], field_eom)
+
 correlations = oqupy.PowerLawSD(alpha=a,
                                 zeta=1,
                                 cutoff=omega_cutoff,
@@ -48,27 +49,29 @@ def H(t, field):
     dt=0.1,
     dkmax=20,
     epsrel=10**(-7))
+
 process_tensor = oqupy.pt_tempo_compute(bath=bath,
                                         start_time=0.0,
                                         end_time=end_time,
                                         parameters=pt_tempo_parameters)
-control = None #oqupy.Control(2)
-dynamics = oqupy.compute_dynamics_single_system_with_field(
-    system=system,
-    process_tensor=process_tensor,
+control = None 
+
+mean_field_dynamics = oqupy.compute_dynamics_with_field(
+    mean_field_system=mean_field_system,
+    process_tensor_list=[process_tensor],
 	initial_field=initial_field,
-    control=control,
+    control_list=[control],
     start_time=0.0,
-    initial_state=initial_state)
-
-print(f"The the final time t = {dynamics.times[-1]:.1f} " \
-      + "the field is {dynamics.fields[-1]:.8g} and the state is:")
-print(dynamics.states[-1])
+    initial_state_list=[initial_state])
+system_dynamics = mean_field_dynamics.system_dynamics[0]
 
+print(f"The final time t = {mean_field_dynamics.times[-1]:.1f} " \
+      + f"the field is {mean_field_dynamics.fields[-1]:.8g} and the state is:")
+print(system_dynamics.states[-1])
 
-t, s_x = dynamics.expectations(oqupy.operators.sigma("x")/2, real=True)
-t, s_z = dynamics.expectations(oqupy.operators.sigma("z")/2, real=True)
-t, fields = dynamics.field_expectations()
+t, s_x = system_dynamics.expectations(oqupy.operators.sigma("x")/2, real=True)
+t, s_z = system_dynamics.expectations(oqupy.operators.sigma("z")/2, real=True)
+t, fields = mean_field_dynamics.field_expectations()
 fig, axes = plt.subplots(2, figsize=(9,10))
 axes[0].plot(t, s_z)
 axes[1].plot(t, np.abs(fields)**2)

examples/mean-field-tempo.py

@@ -25,43 +25,43 @@
 wc=0.0
 kappa=0.1
 end_time=1
-def field_eom(t, state, field):
-    sx_exp = np.matmul(Sx, state).trace().real
+def field_eom(t, states, field):
+    sx_exp = np.matmul(Sx, states[0]).trace().real
     return -(1j*wc+kappa)*field - 1j*gn*sx_exp
-def H(t, field):
-    #print(t, field)
-    return 2.0 * gn * np.abs(field) * Sx
+def H_MF(t, field):
+    return 2.0 * gn * np.real(field) * Sx
 
-system = oqupy.TimeDependentSystemWithField(H,
-        field_eom,
+system = oqupy.TimeDependentSystemWithField(H_MF,
         #gammas = [gam_down, gam_up],
         #lindblad_operators = [oqupy.operators.sigma("-"), tempo.operators.sigma("+")]
         )
+mean_field_system = oqupy.MeanFieldSystem([system], field_eom)
+
 correlations = oqupy.PowerLawSD(alpha=a,
                                 zeta=1,
                                 cutoff=omega_cutoff,
                                 cutoff_type='gaussian',
                                 temperature=temperature)
 bath = oqupy.Bath(oqupy.operators.sigma("z")/2.0, correlations)
 
-
 tempo_parameters = oqupy.TempoParameters(dt=0.1, dkmax=20, epsrel=10**(-7))
 
-tempo_sys = oqupy.TempoWithField(system=system,
-                        bath=bath,
-                        initial_state=initial_state,
+tempo_sys = oqupy.MeanFieldTempo(mean_field_system=mean_field_system,
+                        bath_list=[bath],
+                        initial_state_list=[initial_state],
                         initial_field=initial_field,
                         start_time=0.0,
                         parameters=tempo_parameters)
-dynamics = tempo_sys.compute(end_time=end_time)
+mean_field_dynamics = tempo_sys.compute(end_time=end_time)
+system_dynamics = mean_field_dynamics.system_dynamics[0]
 
-print(f"The the final time t = {dynamics.times[-1]:.1f} " \
-      + f"the field is {dynamics.fields[-1]:.8g} and the state is:")
-print(dynamics.states[-1])
+print(f"The final time t = {mean_field_dynamics.times[-1]:.1f} " \
+      + f"the field is {mean_field_dynamics.fields[-1]:.8g} and the state is:")
+print(system_dynamics.states[-1])
 
-t, s_x = dynamics.expectations(oqupy.operators.sigma("x")/2, real=True)
-t, s_z = dynamics.expectations(oqupy.operators.sigma("z")/2, real=True)
-t, fields = dynamics.field_expectations()
+t, s_x = system_dynamics.expectations(oqupy.operators.sigma("x")/2, real=True)
+t, s_z = system_dynamics.expectations(oqupy.operators.sigma("z")/2, real=True)
+t, fields = mean_field_dynamics.field_expectations()
 fig, axes = plt.subplots(2, figsize=(9,10))
 axes[0].plot(t, s_z)
 axes[1].plot(t, np.abs(fields)**2)