Passing virial back to quantum espresso (fixes #3)

I had to change the place plugin_forces is called in espresso.
This is because in the original implementation plugin_forces
was called *before* stress calculation.
I did some qualitative test, but I am not able to run a proper
thermostated variable cell simulation with QE. Probably will
require some extra validation by QE guys.

I also updated qe version to 5.0.2

patches/qespresso-5.0.1.diff → patches/qespresso-5.0.2.diff

+patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/forces.f90" << \EOF_EOF
+--- ./PW/src/forces.f90.preplumed
++++ ./PW/src/forces.f90
+@@ -247,11 +247,13 @@
+   !
+   IF( textfor ) force(:,:) = force(:,:) + extfor(:,:)
+   !
+   ! ... call void routine for user define/ plugin patches on forces
+   !
+-  CALL plugin_forces()
++  ! plugin should be called after stress has been computed
++  ! Look in pwscf.f90
++  ! CALL plugin_forces()
+   !
+   ! ... write on output the forces
+   !
+   DO na = 1, nat
+      !
+EOF_EOF
 patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_forces.f90" << \EOF_EOF
 --- ./PW/src/plugin_forces.f90.preplumed
 +++ ./PW/src/plugin_forces.f90
-@@ -16,9 +16,42 @@
+@@ -16,9 +16,51 @@
    USE kinds,            ONLY : DP
    USE io_files,         ONLY : outdir
    !
@@ -9,7 +28,7 @@ patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_forces.f90" << \EOF_EOF
    !
 +  USE cell_base,        ONLY : alat, at
 +  USE ions_base,        ONLY : tau, nat,amass
-+  USE force_mod,        ONLY : force
++  USE force_mod,        ONLY : force,sigma
 +  USE control_flags,    ONLY : istep
 +  USE ener,             ONLY : etot 
 +  !
@@ -17,16 +36,22 @@ patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_forces.f90" << \EOF_EOF
    !
 +  INTEGER:: i,j
 +  REAL(DP) :: at_plumed(3,3)
++  REAL(DP) :: virial(3,3)
++  REAL(DP) :: volume
 +  REAL(DP), ALLOCATABLE :: tau_plumed(:,:)
 +  !
-+  ! this is ignored by QE, but required for plumed to work consistently
-+  REAL(DP) :: virial(3,3)
 +  IF(use_plumed) then
 +    IF(ionode)THEN
 +      at_plumed=alat*at;  ! the cell, rescaled properly
 +      allocate(tau_plumed(3,nat))
 +      tau_plumed=alat*tau
-+      virial=0.0
++      volume=+at_plumed(1,1)*at_plumed(2,2)*at_plumed(3,3) &
++             +at_plumed(1,2)*at_plumed(2,3)*at_plumed(3,1) &
++             +at_plumed(1,3)*at_plumed(2,1)*at_plumed(3,2) &
++             -at_plumed(1,1)*at_plumed(3,2)*at_plumed(2,3) &
++             -at_plumed(1,2)*at_plumed(3,3)*at_plumed(2,1) &
++             -at_plumed(1,3)*at_plumed(3,1)*at_plumed(2,2) 
++      virial=-sigma*volume
 +
 +      CALL plumed_f_gcmd("setStep"//char(0),istep)
 +      CALL plumed_f_gcmd("setMasses"//char(0),amass)
@@ -37,9 +62,12 @@ patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_forces.f90" << \EOF_EOF
 +      CALL plumed_f_gcmd("setEnergy"//char(0),etot)
 +      CALL plumed_f_gcmd("calc"//char(0),0)
 +
++      sigma=-virial/volume
++
 +      deallocate(tau_plumed)
 +    ENDIF
 +    CALL mp_bcast(force, ionode_id, intra_image_comm)
++    CALL mp_bcast(sigma, ionode_id, intra_image_comm)
 +  ENDIF
 +  !
    !
@@ -48,7 +76,7 @@ EOF_EOF
 patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_initialization.f90" << \EOF_EOF
 --- ./PW/src/plugin_initialization.f90.preplumed
 +++ ./PW/src/plugin_initialization.f90
-@@ -13,9 +13,49 @@
+@@ -13,9 +13,50 @@
    USE kinds,            ONLY : DP
    USE io_files,         ONLY : tmp_dir
    !
@@ -91,7 +119,8 @@ patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_initialization.f90" <<
 +            call plumed_f_gcmd("setMDEngine"//char(0),"qespresso");
 +            call plumed_f_gcmd("setTimestep"//char(0),dt);
 +            call plumed_f_gcmd("init"//char(0),0);
-+    
++
++
 +        ENDIF
 +     ENDIF
 +  ENDIF
@@ -99,3 +128,22 @@ patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/plugin_initialization.f90" <<
    !
  END SUBROUTINE plugin_initialization
 EOF_EOF
+patch -u -l -b -F 5 --suffix=.preplumed "./PW/src/pwscf.f90" << \EOF_EOF
+--- ./PW/src/pwscf.f90.preplumed
++++ ./PW/src/pwscf.f90
+@@ -126,10 +126,14 @@
+      !
+      ! ... stress calculation
+      !
+      IF ( lstres ) CALL stress ( sigma )
+      !
++     ! ... this is the right place for plugin forces:
++     !
++     CALL plugin_forces()
++     !
+      IF ( lmd .OR. lbfgs ) THEN
+         !
+         if (fix_volume) CALL impose_deviatoric_stress(sigma)
+         !
+         if (fix_area)  CALL  impose_deviatoric_stress_2d(sigma)
+EOF_EOF