Fix multilevel activation for MD (#350)

src/algos/dynamics.f90

@@ -81,14 +81,12 @@ subroutine crest_moleculardynamics(env,tim)
   !>--- init SHAKE? --> we need connectivity info
   if (mddat%shake) then
     calc%calcs(1)%rdwbo = .true.
+    if(.not.calc%calcs(1)%active) calc%calcs(1)%active=.true.
     allocate (grad(3,mol%nat),source=0.0_wp)
     call engrad(mol,calc,energy,grad,io)
     deallocate (grad)
     calc%calcs(1)%rdwbo = .false.
     call move_alloc(calc%calcs(1)%wbo,mddat%shk%wbo)
-    !> moved to within the MD call
-    !call init_shake(mol%nat,mol%at,mol%xyz,mddat%shk,pr)
-    !mddat%nshake = mddat%shk%ncons
   end if
 
   !>--- complete real-time settings to steps

src/calculator/calc_type.f90

@@ -203,6 +203,7 @@ module calc_type
     real(wp),allocatable :: grdtmp(:,:,:)
     real(wp),allocatable :: eweight(:)
     real(wp),allocatable :: weightbackup(:)
+    logical,allocatable  :: activebackup(:)
     real(wp),allocatable :: etmp2(:)
     real(wp),allocatable :: grdtmp2(:,:,:)
     real(wp),allocatable :: eweight2(:)
@@ -640,18 +641,23 @@ subroutine calc_set_active(self,ids)
     integer,intent(in) :: ids(:)
     integer :: i,j,k,l
     if (allocated(self%weightbackup)) deallocate (self%weightbackup)
+    if (allocated(self%activebackup)) deallocate (self%activebackup)
 !>--- on-the-fly multiscale definition
     allocate (self%weightbackup(self%ncalculations),source=1.0_wp)
+    allocate (self%activebackup(self%ncalculations),source=.true.)
     do i = 1,self%ncalculations
 !>--- save backup weights
       self%weightbackup(i) = self%calcs(i)%weight
+      self%activebackup(i) = self%calcs(i)%active
 !>--- set the weight of all unwanted calculations to 0
       if (.not.any(ids(:) .eq. i)) then
         self%calcs(i)%weight = 0.0_wp
         self%calcs(i)%active = .false.
       else
-!>--- and all other to 1
-        self%calcs(i)%weight = 1.0_wp
+!>--- and all other to active
+        if(self%calcs(i)%weight == 0.0_wp)then
+           self%calcs(i)%weight = 1.0_wp
+        endif 
         self%calcs(i)%active = .true.
       end if
     end do
@@ -662,10 +668,11 @@ subroutine calc_set_active_restore(self)
     class(calcdata) :: self
     integer :: i,j,k,l
     if (.not.allocated(self%weightbackup)) return
+    if (.not.allocated(self%activebackup)) return
     do i = 1,self%ncalculations
 !>--- set all to active and restore saved weights
       self%calcs(i)%weight = self%weightbackup(i)
-      self%calcs(i)%active = .true.
+      self%calcs(i)%active = self%activebackup(i)
       self%eweight(i) = self%weightbackup(i)
     end do
     deallocate (self%weightbackup)

src/parsing/confparse2.f90

@@ -109,6 +109,7 @@ subroutine parseinputfile(env,fname)
   call parse_dynamics_data(env,mddat,dict,l1,readstatus)
   if (l1) then
     env%mddat = mddat
+    call env_mddat_specialcases(env)
   end if
 
 !>--- terminate if there were any unrecognized keywords
@@ -212,3 +213,33 @@ subroutine env_calcdat_specialcases(env)
   end if
 
 end subroutine env_calcdat_specialcases
+
+!========================================================================================!
+
+subroutine env_mddat_specialcases(env)
+!****************************************************
+!* Some special treatments are sometimes necessary
+!* depending on a choosen calculation level,
+!* e.g. for on-the-fly multilevel setup
+!* These depend on both the calc data and md data
+!***************************************************
+  use crest_parameters
+  use crest_data
+  use crest_calculator
+  implicit none
+  type(systemdata) :: env
+  integer :: nac,ii,iac
+
+!>--- Check for MD-active only levels
+  if(allocated(env%mddat%active_potentials))then
+    nac = size(env%mddat%active_potentials)
+    do ii=1,nac
+    !>--- deactivate by default (the MD routine will set them to active automatically)  
+     iac = env%mddat%active_potentials(ii)
+     if(iac <= env%calc%ncalculations)then
+       env%calc%calcs(iac)%active = .false.
+     endif
+    enddo
+  endif
+
+end subroutine env_mddat_specialcases