Sandwich Potential

man/xcontrol.7.adoc

@@ -534,6 +534,12 @@ $wall
 *ellipsoid*: auto|'real',auto|'real',auto|'real',all|'list',...::
     set up a ellipsoid wall potential for all or the atoms in 'list'
     with the radii 'real' or an automatical determined sphere radius
+    If 'auto' is chosen for axes, sphere potential is applied, no ellipsoid!!!
+
+*sandwich*: auto|'real',all|'list',...::
+   set up a sandwich wall potential for all or the atoms in 'list'
+   with the radius 'real' or an automatical determined sandwich radius in Bohr
+   Only potential=logfermi ist available. diameter=2*radius+2*4A safety buffer
 
 
 $write

src/constrain_param.f90

@@ -49,6 +49,10 @@
 !>    ellipsoid: auto,<list>
 !>    ellipsoid: <real>,<real>,<real>,all
 !>    ellipsoid: <real>,<real>,<real>,<list>
+!>    sandwich: auto,all
+!>    sandwich: auto,<list>
+!>    sandwich: <real>,all
+!>    sandwich: <real>,<list>
 !> $scan
 !>    ...
 !> $end
@@ -955,8 +959,8 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
    type(TIdentityMap), intent(in) :: idMap
    real(wp),intent(in) :: xyz(3,nat)
 
-   integer  :: idum
-   real(wp) :: ddum,darray(3)
+   integer  :: idum,i
+   real(wp) :: ddum,darray(3),min_z,max_z
    logical  :: ldum
    integer  :: list(nat),nlist
    integer  :: tlist(nat),ntlist
@@ -985,7 +989,7 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
    case default ! ignore, don't even think about raising them
    case('sphere')
       if (narg.lt.2) then
-         call env%warning("Not enough arguments to set up a spherical wall",source)
+         call env%error("Not enough arguments to set up a spherical wall",source)
          return
       endif
    !  part 1: get the sphere radius
@@ -1007,11 +1011,11 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
          do iarg = 2, narg
             if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then
                if (nlist+ntlist.gt.nat) then
-                  call env%warning("Too many atoms in list for spherical wall.",source)
+                  call env%error("Too many atoms in list for spherical wall.",source)
                   return ! something went wrong
                endif
                if (maxval(tlist(:ntlist)).gt.nat) then
-                  call env%warning("Attempted to wall in a non-existing atom",source)
+                  call env%error("Attempted to wall in a non-existing atom",source)
                   cycle ! skip crappy input
                endif
                list(nlist+1:nlist+ntlist) = tlist
@@ -1023,12 +1027,12 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
          enddo
          call set_sphere_radius(radius,center,nlist,list)
       endif
-      write(env%unit,'("spherical wallpotenial with radius",'//&
+      write(env%unit,'("spherical wallpotential with radius",'//&
          '1x,f12.7,1x,"Å")') radius*autoaa
 
    case('ellipsoid')
       if (narg.lt.4) then
-         call env%warning("Not enough arguments to set up an ellipsoidal wall",source)
+         call env%error("Not enough arguments to set up an ellipsoidal wall",source)
          return
       endif
    !  part 1: get ellipsoid axis
@@ -1057,11 +1061,11 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
          do iarg = 4, narg
             if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then
                if (nlist+ntlist.gt.nat) then
-                  call env%warning("Too many atoms in list for spherical wall.",source)
+                  call env%error("Too many atoms in list for spherical wall.",source)
                   return ! something went wrong
                endif
                if (maxval(tlist(:ntlist)).gt.nat) then
-                  call env%warning("Attempted to wall in a non-existing atom",source)
+                  call env%error("Attempted to wall in a non-existing atom",source)
                   cycle ! skip crappy input
                endif
                list(nlist+1:nlist+ntlist) = tlist
@@ -1076,6 +1080,70 @@ subroutine set_wall(env,key,val,nat,at,idMap,xyz)
       write(env%unit,'("ellipsoidal wallpotenial with radii",'//&
          '3(1x,f12.7,1x,"Å"))') radii*autoaa
 
+   case('sandwich')
+      if (narg.lt.2) then
+         call env%error("Not enough arguments to set up sandwich walls",source)
+         return
+      endif
+   !  part 1: get the sandwich distance
+      wpot%sandwich = .true.
+      set%do_cma_trafo = .true.
+      center = 0.0_wp
+!      number_walls=1
+      call get_sphere_radius(nat,at,xyz,center,radius,do_trafo=.true.)
+      if (trim(argv(1)).eq.'auto') then
+            radius=(maxval(xyz(3,:))-minval(xyz(3,:)))/2.0_wp ! need to set $cma in xcontrol, not done automatically
+            wpot(1)%radius=radius
+      else
+            if (getValue(env,trim(argv(1)),ddum)) then
+                radius = ddum !in Bohr!!!
+                wpot(1)%radius=radius
+            else
+                call env%error("Undefined arguments for sandwich: ... in your xcontrol file!",source)
+                return ! something went wrong
+            endif
+      endif
+
+   !  part 2: get atoms 
+      if (trim(argv(2)).eq.'all') then
+         call set_sphere_radius(radius,center)
+      else
+         do iarg = 2, narg
+            if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then
+               if (nlist+ntlist.gt.nat) then
+                  call env%error("Too many atoms in list for spherical wall.",source)
+                  return ! something went wrong
+               endif
+               if (maxval(tlist(:ntlist)).gt.nat) then
+                  call env%error("Attempted to wall in a non-existing atom.",source)
+                  cycle ! skip crappy input
+               endif
+               list(nlist+1:nlist+ntlist) = tlist
+               nlist = nlist + ntlist
+
+               !get auto sandwich distance for list of atoms
+               max_z = 0.0_wp
+               min_z = 0.0_wp
+               do i = 1, nat
+                  if (any(list == i)) then
+                     max_z = max(max_z, xyz(3,i))
+                     min_z = min(min_z, xyz(3,i))
+                  end if
+               end do
+               radius=(max_z - min_z)/2.0_wp
+               wpot(1)%radius=radius
+
+            else
+               ! warning already generated by get_list_value
+               return ! something went wrong
+            endif
+         enddo
+         call set_sphere_radius(radius,center,nlist,list)
+      endif
+
+      write(env%unit,'("sandwich wallpotential with radius in A (diameter=2*radius+2*4A safety buffer) ",'//&
+         '1x,f12.7,1x,"Å")') radius*autoaa
+
    end select
 
 end subroutine set_wall
@@ -1513,6 +1581,7 @@ subroutine set_legacy(env,key,val,nat,at,idMap,xyz)
 !   case('scan')
    case('ellips'); call set_wall(env,'ellipsoid',val,nat,at,idMap,xyz)
    case('sphere'); call set_wall(env,'sphere',val,nat,at,idMap,xyz)
+   case('sandwich'); call set_wall(env,'sandwich',val,nat,at,idMap,xyz)
    case('fix'); call set_fix(env,'atoms',val,nat,at,idMap,xyz)
    case('atomlist+'); call set_metadyn(env,'atoms',val,nat,at,idMap,xyz)
    end select

src/gfnff/calculator.f90

@@ -256,7 +256,6 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
    ! ---------------------------------------!
    ! post processing of gradient and energy !
    !----------------------------------------!
-
    ! various external potentials !
    call constrain_pot(potset,mol%n,mol%at,mol%xyz,gradient,efix)
    call constrpot   (mol%n,mol%at,mol%xyz,gradient,efix)

src/prog/main.F90

@@ -377,7 +377,10 @@ subroutine xtbMain(env, argParser)
          call init_metadyn(mol%n, metaset%maxsave)
       end if
       call load_rmsdbias(rmsdset, mol%n, mol%at, mol%xyz)
-
+      ! ------------------------------------------------------------------------
+      !> CONSTRAINTS & SCANS
+      !> now we are at a point that we can check for requested constraints
+      call read_userdata(xcontrol, env, mol)
       ! ------------------------------------------------------------------------
       !> get some memory
       allocate (cn(mol%n), sat(mol%n), g(3, mol%n), source=0.0_wp)
@@ -420,10 +423,6 @@ subroutine xtbMain(env, argParser)
          mol%info%two_dimensional = .false.
       end if
 
-      ! ------------------------------------------------------------------------
-      !> CONSTRAINTS & SCANS
-      !> now we are at a point that we can check for requested constraints
-      call read_userdata(xcontrol, env, mol)
 
       !> initialize metadynamics
       call load_metadynamic(metaset, mol%n, mol%at, mol%xyz)