PINT: Multiple-time striding in imaginary time to couple PIMD solute …

…to fewer PIMC solvent beads
cp2k · Sep 3, 2019 · 2e6ac99 · 2e6ac99
1 parent ae9949d
commit 2e6ac99
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 6 deletions.
diff --git a/src/motion/helium_interactions.F b/src/motion/helium_interactions.F
@@ -312,7 +312,7 @@ SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
       CHARACTER(len=*), PARAMETER :: routineN = 'helium_bead_solute_e_f', &
          routineP = moduleN//':'//routineN
 
-      INTEGER                                            :: hbeads, hi, qi
+      INTEGER                                            :: hbeads, hi, qi, stride
       REAL(KIND=dp), DIMENSION(3)                        :: helium_r
       REAL(KIND=dp), DIMENSION(:), POINTER               :: my_force
 
@@ -361,6 +361,15 @@ SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
 
       END SELECT
 
+      ! Account for Imaginary time striding in forces:
+      IF (PRESENT(force)) THEN
+         stride = 1
+         IF (hbeads < pint_env%p) THEN
+            stride = pint_env%p/hbeads
+         END IF
+         force = force*REAL(stride,dp)
+      END IF
+
       RETURN
    END SUBROUTINE helium_bead_solute_e_f
 

diff --git a/src/motion/helium_methods.F b/src/motion/helium_methods.F
@@ -633,11 +633,23 @@ SUBROUTINE helium_create(helium_env, input, solute)
             ! fill in the solute-related data structures
             helium_env(k)%helium%e_corr = 0.0_dp
             IF (helium_env(k)%helium%solute_present) THEN
-               helium_env(k)%helium%bead_ratio = helium_env(k)%helium%beads/helium_env(k)%helium%solute_beads
-
-               ! check if bead numbers are commensurate:
-               i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%solute_beads-helium_env(k)%helium%beads
-               CPASSERT(i == 0)
+               IF (helium_env(k)%helium%solute_beads > helium_env(k)%helium%beads) THEN
+                  ! Imaginary time striding for solute:
+                  helium_env(k)%helium%bead_ratio = helium_env(k)%helium%solute_beads/&
+                     helium_env(k)%helium%beads
+                  ! check if bead numbers are commensurate:
+                  i = helium_env(k)%helium%bead_ratio*helium_env(k)%helium%beads-helium_env(k)%helium%solute_beads
+                  IF (i /= 0) THEN
+                     msg_str = "Adjust number of solute beads to multiple of solvent beads."
+                     CPABORT(msg_str)
+                  END IF
+                  msg_str = "Using multiple-time stepping in imaginary time for solute to couple "// &
+                     "to fewer solvent beads. Warning: Avoid too large coupling factors."
+                  CPWARN(msg_str)
+               ELSE IF  (helium_env(k)%helium%solute_beads < helium_env(k)%helium%beads)
+                  msg_str = "Solute beads < solvent beads: Not supported."
+                  CPABORT(msg_str)
+               END IF
 !TODO       Adjust helium bead number if not comm. and if coords not given expl.
 
                ! check if tau, temperature and bead number are consistent: