1! This file is part of xtb.
2!
3! Copyright (C) 2017-2020 Stefan Grimme
4!
5! xtb is free software: you can redistribute it and/or modify it under
6! the terms of the GNU Lesser General Public License as published by
7! the Free Software Foundation, either version 3 of the License, or
8! (at your option) any later version.
9!
10! xtb is distributed in the hope that it will be useful,
11! but WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13! GNU Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public License
16! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
17
18!! ========================================================================
19!  * WELCOME TO THE   C O N S T R A I N S   &   S C A N S   MODULE IN XTB *
20!! ------------------------------------------------------------------------
21!  The syntax is:
22!> $constrain
23!> ## this part is read by setparam.f90
24!>    force constant=<real>
25!>    all bonds=<bool>
26!>    all angles=<bool>
27!>    all torsions=<bool>
28!> ## this part is read in this module
29!>    distance: <i>,<j>,        auto
30!>    distance: <i>,<j>,        <real>
31!>    angle:    <i>,<j>,<k>,    auto
32!>    angle:    <i>,<j>,<k>,    <real>
33!>    dihedral: <i>,<j>,<k>,<l>,auto
34!>    dihedral: <i>,<j>,<k>,<l>,<real>
35!>    center:   <int>,<real>
36!> $fix
37!>    force constant=<real>
38!>    spring exponent=<int>
39!>    atoms: <list>
40!> $split
41!>    fragment1: <list>
42!>    fragment2: <list>
43!> $wall
44!>    sphere: auto,all
45!>    sphere: auto,<list>
46!>    sphere: <real>,all
47!>    sphere: <real>,<list>
48!>    ellipsoid: auto,all
49!>    ellipsoid: auto,<list>
50!>    ellipsoid: <real>,<real>,<real>,all
51!>    ellipsoid: <real>,<real>,<real>,<list>
52!> $scan
53!>    ...
54!> $end
55!! ========================================================================
56module xtb_constrain_param
57   use xtb_mctc_accuracy, only : wp
58   use xtb_mctc_strings, only : parse
59   use xtb_readin, only : getline => strip_line,getValue,getListValue
60   use xtb_setparam, only : verbose
61   use xtb_type_environment, only : TEnvironment
62   use xtb_type_identitymap, only : TIdentityMap, init
63   use xtb_type_molecule, only : TMolecule
64
65   implicit none
66
67   private :: getline,getValue, handlerInterface
68
69   character,private,parameter :: flag = '$'
70   character,private,parameter :: colon = ':'
71   character,private,parameter :: space = ' '
72   character,private,parameter :: equal = '='
73   character,private,parameter :: hash = '#'
74   character,private,parameter :: comma = ','
75   character,private,parameter :: semicolon = ';'
76   character,private,parameter :: dot = '.'
77   character(len=*),private,parameter :: flag_end = flag//'end'
78
79!  Using allocatable arrays of dynamic length strings is only possible
80!  with a lot of hacks, so we use good'ol fixed size stack arrays.
81!  Let's choose something different from 42 that is not dividable by 10... ;)
82!  Happy debugging!
83   integer,private,parameter :: p_str_length = 48
84   integer,private,parameter :: p_arg_length = 24
85
86   public
87
88   abstract interface
89      subroutine handlerInterface(env,key,val,nat,at,idMap,xyz)
90         import :: wp, TEnvironment, TIdentityMap
91         type(TEnvironment), intent(inout) :: env
92         character(len=*),intent(in) :: key
93         character(len=*),intent(in) :: val
94         integer, intent(in) :: nat
95         type(TIdentityMap), intent(in) :: idMap
96         integer, intent(in) :: at(nat)
97         real(wp),intent(in) :: xyz(3,nat)
98      end subroutine handlerInterface
99   end interface
100
101contains
102
103subroutine read_userdata(fname,env,mol)
104   use xtb_readin, only : find_new_name
105   use xtb_scanparam
106   implicit none
107   character(len=*), parameter :: source = 'userdata_read'
108   type(TEnvironment), intent(inout) :: env
109   type(TMolecule), intent(inout) :: mol
110   character(len=*),intent(in)  :: fname
111   character(len=:),allocatable :: line
112   character(len=:),allocatable :: key
113   character(len=:),allocatable :: val
114   character(len=:),allocatable :: newname
115   type(TIdentityMap) :: idMap
116   integer :: i
117   integer :: id
118   integer :: ic
119   integer :: ie
120   integer :: err
121   logical :: exist
122
123   if (verbose) then
124      write(env%unit,'(72("$"))')
125      write(env%unit,'(1x,"CONSTRAINTS & SCANS: DEBUG SECTION")')
126      write(env%unit,'(72("$"))')
127   endif
128
129   call open_file(id,fname,'r')
130   if (id.eq.-1) then
131      call env%warning("could not find '"//fname//"'",source)
132      return
133   endif
134   rewind(id) ! not sure if this is necessary
135
136   call init(idMap, mol)
137
138   call idMap%writeInfo(env%unit)
139
140!  read first line before the readloop starts, I have to do this
141!  to avoid using backspace on id (dammit Turbomole format)
142   call getline(id,line,err)
143   readflags: do
144      !  check if there is a $ in the *first* column
145      if (index(line,flag).eq.1) then
146         select case(line(2:))
147         case('fix'      )
148            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
149            call rdblock(env,set_fix,    line,id,mol%n,mol%at,idMap,mol%xyz,err)
150         case('split'    )
151            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
152            call rdblock(env,set_split,  line,id,mol%n,mol%at,idMap,mol%xyz,err)
153         case('constrain')
154            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
155            if (allocated(potset%xyz)) then
156               call rdblock(env,set_constr, line,id,mol%n,mol%at,idMap,potset%xyz,err)
157            else
158               call rdblock(env,set_constr, line,id,mol%n,mol%at,idMap,mol%xyz,err)
159            endif
160         case('scan'     )
161            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
162            call rdblock(env,set_scan,   line,id,mol%n,mol%at,idMap,mol%xyz,err)
163         case('wall'     )
164            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
165            call rdblock(env,set_wall,   line,id,mol%n,mol%at,idMap,mol%xyz,err)
166         case('metadyn'  )
167            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
168            call rdblock(env,set_metadyn,line,id,mol%n,mol%at,idMap,mol%xyz,err)
169         case('hess'     )
170            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
171            call rdblock(env,set_hess,   line,id,mol%n,mol%at,idMap,mol%xyz,err)
172         case('path'     )
173            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
174            call rdblock(env,set_path,   line,id,mol%n,mol%at,idMap,mol%xyz,err)
175         case('reactor'  )
176            if (verbose) write(env%unit,'(">",1x,a)') line(2:)
177            call rdblock(env,set_reactor,line,id,mol%n,mol%at,idMap,mol%xyz,err)
178         case('set'      ); call rdsetbl(env,set_legacy,line,id,mol%n,mol%at,idMap,mol%xyz,err)
179         case default ! unknown keyword -> ignore, we don't raise them
180            call getline(id,line,err)
181         end select
182      else ! not a keyword -> ignore
183         call getline(id,line,err)
184      endif
185   !  check for end of file, which I will tolerate as alternative to $end
186      if (is_iostat_end(err)) exit readflags
187!     if (index(line,flag_end).ne.0) exit readflags ! compatibility reasons
188   enddo readflags
189
190   if (verbose) write(env%unit,'(72("$"))')
191   call close_file(id)
192end subroutine read_userdata
193
194subroutine rdsetbl(env,handler,line,id,nat,at,idMap,xyz,err)
195   implicit none
196   character(len=*), parameter :: source = 'userdata_rdsetbl'
197   type(TEnvironment), intent(inout) :: env
198   integer,intent(in) :: id
199   procedure(handlerInterface) :: handler
200   integer, intent(in) :: nat
201   integer, intent(in) :: at(nat)
202   type(TIdentityMap), intent(in) :: idMap
203   real(wp),intent(in) :: xyz(3,nat)
204   integer,intent(out) :: err
205   character(len=:),allocatable,intent(out) :: line
206   character(len=:),allocatable :: key
207   character(len=:),allocatable :: val
208   integer :: ie
209   logical :: exitRun
210   do
211      call getline(id,line,err)
212      if (is_iostat_end(err)) exit
213      if (index(line,flag).ne.0) exit
214      if (verbose) write(env%unit,'("->",1x,a)') line
215
216      ! find the first colon
217      ie = index(line,space)
218      if ((line.eq.'').or.(ie.eq.0)) cycle
219      key = trim(line(:ie-1))
220      val = trim(adjustl(line(ie+1:)))
221      call handler(env,key,val,nat,at,idMap,xyz)
222      call env%check(exitRun)
223      if (exitRun) then
224         call env%error("handler could not process input", source)
225         return
226      end if
227   enddo
228
229end subroutine rdsetbl
230
231subroutine rdblock(env,handler,line,id,nat,at,idMap,xyz,err)
232   implicit none
233   character(len=*), parameter :: source = 'userdata_rdblock'
234   type(TEnvironment), intent(inout) :: env
235   integer,intent(in) :: id
236   procedure(handlerInterface) :: handler
237   integer, intent(in) :: nat
238   integer, intent(in) :: at(nat)
239   type(TIdentityMap), intent(in) :: idMap
240   real(wp),intent(in) :: xyz(3,nat)
241   integer,intent(out) :: err
242   character(len=:),allocatable,intent(out) :: line
243   character(len=:),allocatable :: key
244   character(len=:),allocatable :: val
245   integer :: ie
246   logical :: exitRun
247   do
248      call getline(id,line,err)
249      if (is_iostat_end(err)) exit
250      if (index(line,flag).ne.0) exit
251      if (verbose) write(env%unit,'("->",1x,a)') line
252
253      ! find the first colon
254      ie = index(line,colon)
255      if ((line.eq.'').or.(ie.eq.0)) cycle
256      key = trim(line(:ie-1))
257      val = trim(adjustl(line(ie+1:)))
258      call handler(env,key,val,nat,at,idMap,xyz)
259      call env%check(exitRun)
260      if (exitRun) then
261         call env%error("handler could not process input", source)
262         return
263      end if
264   enddo
265
266end subroutine rdblock
267
268subroutine set_fix(env,key,val,nat,at,idMap,xyz)
269   use xtb_type_atomlist
270   use xtb_fixparam
271   use xtb_setparam
272   implicit none
273   character(len=*), parameter :: source = 'userdata_fix'
274   type(TEnvironment), intent(inout) :: env
275   character(len=*),intent(in) :: key
276   character(len=*),intent(in) :: val
277   integer, intent(in) :: nat
278   integer, intent(in) :: at(nat)
279   type(TIdentityMap), intent(in) :: idMap
280   real(wp),intent(in) :: xyz(3,nat)
281
282   type(TAtomList) :: atl
283
284   integer  :: i
285   integer  :: iat
286   integer  :: idum
287   integer  :: nlist
288   integer, allocatable :: list(:)
289   real(wp) :: ddum
290   logical  :: ldum
291
292   integer  :: narg
293   character(len=p_str_length),dimension(p_arg_length) :: argv
294
295   call atl%resize(nat)
296
297   call parse(val,comma,argv,narg)
298!  some debug xtb_printout
299   if (verbose) then
300      do idum = 1, narg
301         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
302      enddo
303   endif
304   select case(key)
305   case default ! ignore, don't even think about raising them
306   case('elements')
307      call atl%new
308      do idum = 1, narg
309         ! get element by symbol
310         if (idMap%has(argv(idum))) then
311            call idMap%get(list, argv(idum))
312            if (allocated(list)) then
313               call atl%add(list)
314            else
315               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
316               cycle
317            end if
318         else
319            ldum = getValue(env,trim(argv(idum)),iat)
320            if (.not.ldum) cycle ! skip garbage input
321            ! check for unreasonable input
322            if (iat > 0) then
323               call atl%add(at.eq.iat)
324            else
325               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
326               cycle
327            endif
328         endif
329      enddo
330      if (fixset%n > 0) call atl%add(fixset%atoms(:fixset%n))
331      call atl%to_list(list)
332      fixset%atoms = list
333      fixset%n = size(list)
334   case('atoms')
335      call atl%new(val)
336      if (atl%get_error()) then
337         call env%warning('something is wrong in the fixing list',source)
338         return
339      endif
340      if (fixset%n > 0) call atl%add(fixset%atoms(:fixset%n))
341      call atl%to_list(list)
342      fixset%atoms = list
343      fixset%n = size(list)
344   case('freeze')
345      call atl%new(val)
346      if (atl%get_error()) then
347         call env%warning('something is wrong in the freezing list',source)
348         return
349      endif
350      if (freezeset%n > 0) call atl%add(freezeset%atoms(:freezeset%n))
351      call atl%to_list(list)
352      freezeset%atoms = list
353      freezeset%n = size(list)
354   case('shake')
355      allocate(list(nat*(nat+1)/2), source=0)
356      if (mod(narg,2).ne.0) then
357         call env%warning("could not read input for user defined shake!",source)
358         return
359      endif
360         if (narg+shakeset%n > nat*(nat+1)/2) then
361         call env%warning("too many SHAKE constraints!",source)
362         return
363      endif
364      if (.not.shake_md) shake_md = .true.
365      do idum = 1, narg
366         if (getValue(env,trim(argv(idum)),iat)) then
367            if (iat.gt.nat) then
368               call env%warning('Attempted constrain atom not present in molecule.',source)
369               cycle
370            endif
371            shakeset%n = shakeset%n+1
372            shakeset%atoms(shakeset%n) = iat
373         else
374            call env%warning("Something went wrong in set_fix_ 'shake'.",source)
375            return ! you screwed it, let's get out of here
376         endif
377      enddo
378   end select
379
380end subroutine set_fix
381
382subroutine set_constr(env,key,val,nat,at,idMap,xyz)
383   use xtb_mctc_constants
384   use xtb_mctc_convert
385   use xtb_type_atomlist
386   use xtb_scanparam
387   use xtb_splitparam
388   implicit none
389   character(len=*), parameter :: source = 'userdata_constr'
390   type(TEnvironment), intent(inout) :: env
391   character(len=*),intent(in) :: key
392   character(len=*),intent(in) :: val
393   integer, intent(in) :: nat
394   integer, intent(in) :: at(nat)
395   type(TIdentityMap), intent(in) :: idMap
396   real(wp),intent(in) :: xyz(3,nat)
397
398   type(TAtomList) :: atl
399
400   integer  :: iat
401   integer  :: ioffset
402   integer  :: idum
403   real(wp) :: ddum
404   integer  :: nlist
405   integer, allocatable :: list(:)
406   logical  :: ldum
407   integer  :: i,j,k,l
408   real(wp) :: phi,dist,ra(3),rb(3)
409   real(wp),external :: valijkl
410
411   integer  :: narg
412   character(len=p_str_length),dimension(p_arg_length) :: argv
413
414   call atl%resize(nat)
415
416   call parse(val,comma,argv,narg)
417   if (verbose) then
418      do idum = 1, narg
419         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
420      enddo
421   endif
422   select case(key)
423   case default ! ignore, don't even think about raising them
424
425   case('elements')
426      call atl%new
427      do idum = 1, narg
428         ! get element by symbol
429         if (idMap%has(argv(idum))) then
430            call idMap%get(list, argv(idum))
431            if (allocated(list)) then
432               call atl%add(list)
433            else
434               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
435               cycle
436            end if
437         else
438            ldum = getValue(env,trim(argv(idum)),iat)
439            if (.not.ldum) cycle ! skip garbage input
440            ! check for unreasonable input
441            if (iat > 0) then
442               call atl%add(at.eq.iat)
443            else
444               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
445               cycle
446            endif
447         endif
448      enddo
449      if (potset%pos%n > 0) call atl%add(potset%pos%atoms(:potset%pos%n))
450      call atl%to_list(list)
451      potset%pos%atoms = list
452      potset%pos%n = size(list)
453   case('atoms')
454      call atl%new(val)
455      if (atl%get_error()) then
456         call env%warning('something is wrong in the fixing list',source)
457         return
458      endif
459      if (potset%pos%n > 0) call atl%add(potset%pos%atoms(:potset%pos%n))
460      call atl%to_list(list)
461      potset%pos%atoms = list
462      potset%pos%n = size(list)
463
464   case('DISTANCE')
465      if (narg.lt.3 .or. narg.gt.4) then
466         call env%error('not enough arguments to constrain a distance',source)
467         return
468      endif
469      ioffset = 2*potset%dist%n
470      potset%dist%n = potset%dist%n+1
471   !  part 1: get the constrained atoms
472      do i = 1, 2
473         if (getValue(env,trim(argv(i)),idum)) then
474            potset%dist%atoms(ioffset+i) = idum
475         else
476            call env%warning("Something went wrong in set_constr_ 'distance'. (1)",source)
477            potset%dist%n = potset%dist%n-1 ! remove invalid contrain
478            return
479         endif
480      enddo
481   !  part 2: get the distance between those atoms
482      i = potset%dist%atoms(ioffset+1)
483      j = potset%dist%atoms(ioffset+2)
484      dist = sqrt(sum((xyz(:,i)-xyz(:,j))**2))
485      if (trim(argv(3)).eq.'auto') then
486         potset%dist%val(potset%dist%n) = dist
487      else
488         if (getValue(env,trim(argv(3)),ddum)) then
489            potset%dist%val(potset%dist%n) = ddum * aatoau
490         else
491            call env%warning("Something went wrong in set_constr_ 'distance'. (2)",source)
492            potset%dist%n = potset%dist%n-1 ! remove invalid contrain
493            return
494         endif
495      endif
496      if (narg.eq.4) then
497         if (getValue(env,trim(argv(4)),idum)) then
498            if (idum < 2 .or. mod(idum, 2).ne.0) then
499               call env%warning("Invalid spring exponent given", source)
500               potset%dist%n = potset%dist%n-1 ! remove invalid contrain
501               return
502            end if
503            potset%dist%expo(potset%dist%n) = real(idum, wp)
504         else
505            call env%warning("Something went wrong in set_constr_ 'distance'. (3)",source)
506            potset%dist%n = potset%dist%n-1 ! remove invalid contrain
507            return
508         end if
509         write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//&
510            '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å",1x,"with expo",1x,i0)') &
511            i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa, &
512            nint(potset%dist%expo(potset%dist%n))
513      else
514         write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//&
515            '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å")') &
516            i,j, potset%dist%val(potset%dist%n)*autoaa, dist*autoaa
517      end if
518
519   case('ANGLE')
520      if (narg.ne.4) then
521         call env%error('not enough arguments to constrain an angle',source)
522         return
523      endif
524      ioffset = 3*potset%angle%n
525      potset%angle%n = potset%angle%n+1
526   !  part 1: get the constrained atoms
527      do i = 1, narg-1
528         if (getValue(env,trim(argv(i)),idum)) then
529            potset%angle%atoms(ioffset+i) = idum
530         else
531            call env%warning("Something went wrong in set_constr_ 'angle'. (1)",source)
532            potset%angle%n = potset%angle%n-1 ! remove invalid contrain
533            return
534         endif
535      enddo
536   !  part 2: get the angle between the constrained atoms
537      i = potset%angle%atoms(ioffset+1)
538      j = potset%angle%atoms(ioffset+2)
539      k = potset%angle%atoms(ioffset+3)
540      call bangl(xyz,i,j,k,phi)
541      if (trim(argv(narg)).eq.'auto') then
542         potset%angle%val(potset%angle%n) = phi
543      else
544         if (getValue(env,trim(argv(narg)),ddum)) then
545            potset%angle%val(potset%angle%n) = pi/180.0_wp * ddum
546         else
547            call env%warning("Something went wrong in set_constr_ 'angle'. (2)",source)
548            potset%angle%n = potset%angle%n-1 ! remove invalid contrain
549            return
550         endif
551      endif
552
553   case('DIHEDRAL')
554      if (narg.ne.5) then
555         call env%error('not enough arguments to constrain a dihedral',source)
556         return
557      endif
558      ioffset = 4*potset%dihedral%n
559      potset%dihedral%n = potset%dihedral%n+1
560      if (nconstr.gt.maxconstr) then ! double check this
561         call env%error('This should never happen! Let somebody check set_constr',source)
562         return
563      end if
564   !  part 1: get the constrained atoms
565      do i = 1, narg-1
566         if (getValue(env,trim(argv(i)),idum)) then
567            potset%dihedral%atoms(ioffset+i) = idum
568         else
569            call env%warning("Something went wrong in set_constr_ 'dihedral'. (1)",source)
570            potset%dihedral%n = potset%dihedral%n-1 ! remove invalid contrain
571            return
572         endif
573      enddo
574   !  part 2: get the angle between the constrained atoms
575      i = potset%dihedral%atoms(ioffset+1)
576      j = potset%dihedral%atoms(ioffset+2)
577      k = potset%dihedral%atoms(ioffset+3)
578      l = potset%dihedral%atoms(ioffset+4)
579      phi=valijkl(nat,xyz,i,j,k,l)
580      if (trim(argv(narg)).eq.'auto') then
581         potset%dihedral%val(potset%dihedral%n) = phi
582      else
583         if (getValue(env,trim(argv(narg)),ddum)) then
584            potset%dihedral%val(potset%dihedral%n) = pi/180.0_wp * ddum
585         else
586            call env%warning("Something went wrong in set_constr_ 'dihedral'. (2)",source)
587            potset%dihedral%n = potset%dihedral%n-1 ! remove invalid contrain
588            return
589         endif
590      endif
591
592   case('distance')
593      if (narg.ne.3) then
594         call env%error('not enough arguments to constrain a distance',source)
595         return
596      endif
597      nconstr = nconstr+1
598      if (nconstr.gt.maxconstr) then ! double check this
599         call env%error('This should never happen! Let somebody check set_constr',source)
600         return
601      end if
602   !  part 1: get the constrained atoms
603      do i = 1, narg-1
604         if (getValue(env,trim(argv(i)),idum)) then
605            atconstr(i,nconstr) = idum
606         else
607            call env%warning("Something went wrong in set_constr_ 'distance'. (1)",source)
608            nconstr = nconstr-1 ! remove invalid contrain
609            return
610         endif
611      enddo
612   !  part 2: get the distance between those atoms
613      i = atconstr(1,nconstr)
614      j = atconstr(2,nconstr)
615      dist = sqrt(sum((xyz(:,i)-xyz(:,j))**2))
616      if (trim(argv(narg)).eq.'auto') then
617         valconstr(nconstr) = dist
618      else
619         if (getValue(env,trim(argv(narg)),ddum)) then
620            valconstr(nconstr) = ddum * aatoau
621         else
622            call env%warning("Something went wrong in set_constr_ 'distance'. (2)",source)
623            nconstr = nconstr-1 ! remove invalid contrain
624            return
625         endif
626      endif
627      write(env%unit,'("constraining bond",2(1x,i0),1x,"to",'//&
628         '1x,f12.7,1x,"Å, actual value:",1x,f12.7,1x,"Å")') &
629         i,j, valconstr(nconstr)*autoaa, dist*autoaa
630
631   case('angle')
632      if (narg.ne.4) then
633         call env%error('not enough arguments to constrain an angle',source)
634         return
635      endif
636      nconstr = nconstr+1
637      if (nconstr.gt.maxconstr) then ! double check this
638         call env%error('This should never happen! Let somebody check set_constr',source)
639         return
640      end if
641   !  part 1: get the constrained atoms
642      do i = 1, narg-1
643         if (getValue(env,trim(argv(i)),idum)) then
644            atconstr(i,nconstr) = idum
645         else
646            call env%warning("Something went wrong in set_constr_ 'angle'. (1)",source)
647            nconstr = nconstr-1 ! remove invalid contrain
648            return
649         endif
650      enddo
651   !  part 2: get the angle between the constrained atoms
652      i = atconstr(1,nconstr)
653      j = atconstr(2,nconstr)
654      k = atconstr(3,nconstr)
655      call bangl(xyz,i,j,k,phi)
656      if (trim(argv(narg)).eq.'auto') then
657         valconstr(nconstr) = phi
658      else
659         if (getValue(env,trim(argv(narg)),ddum)) then
660            valconstr(nconstr) = pi/180.0_wp * ddum
661         else
662            call env%warning("Something went wrong in set_constr_ 'angle'. (2)",source)
663            nconstr = nconstr-1 ! remove invalid contrain
664            return
665         endif
666      endif
667      write(env%unit,'("constraining angle",3(1x,i0),1x,"to",'//&
668         '1x,f12.7,"°, actual value:",1x,f12.7,"°")') &
669         i,j,k,180.0_wp/pi * valconstr(nconstr),180.0_wp/pi * phi
670
671   case('dihedral')
672      if (narg.ne.5) then
673         call env%error('not enough arguments to constrain a dihedral',source)
674         return
675      endif
676      nconstr = nconstr+1
677      if (nconstr.gt.maxconstr) then ! double check this
678         call env%error('This should never happen! Let somebody check set_constr',source)
679         return
680      end if
681   !  part 1: get the constrained atoms
682      do i = 1, narg-1
683         if (getValue(env,trim(argv(i)),idum)) then
684            atconstr(i,nconstr) = idum
685         else
686            call env%warning("Something went wrong in set_constr_ 'dihedral'. (1)",source)
687            nconstr = nconstr-1 ! remove invalid contrain
688            return
689         endif
690      enddo
691   !  part 2: get the angle between the constrained atoms
692      i = atconstr(1,nconstr)
693      j = atconstr(2,nconstr)
694      k = atconstr(3,nconstr)
695      l = atconstr(4,nconstr)
696      phi=valijkl(nat,xyz,i,j,k,l)
697      if (trim(argv(narg)).eq.'auto') then
698         valconstr(nconstr) = phi
699      else
700         if (getValue(env,trim(argv(narg)),ddum)) then
701            valconstr(nconstr) = pi/180.0_wp * ddum
702         else
703            call env%warning("Something went wrong in set_constr_ 'dihedral'. (2)",source)
704            nconstr = nconstr-1 ! remove invalid contrain
705            return
706         endif
707      endif
708      write(env%unit,'("constraining angle",4(1x,i0),1x,"to",'//&
709         '1x,f12.7,"°, actual value:",1x,f12.7,"°")') &
710         i,j,k,l,180.0_wp/pi * valconstr(nconstr),180.0_wp/pi * phi
711
712   case('center')
713      ! copied from molbld.f without modification (except for error handling)
714      if (narg.ne.2) then
715         call env%error('not enough argument to constrain center',source)
716         return
717      endif
718      nconstr = nconstr+1
719      atconstr(1,nconstr) = -2
720      if (getValue(env,trim(argv(2)),ddum)) then
721         valconstr(nconstr) = ddum
722      else
723         call env%warning("Something went wrong in set_constr_ 'center'. (1)",source)
724         nconstr = nconstr-1
725         return
726      endif
727      if (getValue(env,trim(argv(1)),idum)) then
728         iatf1 = idum
729      else
730         call env%warning("Something went wrong in set_constr_ 'center'. (2)",source)
731         nconstr = nconstr-1
732         return
733      endif
734      massf1 = 0.0_wp
735      do i = 1, iatf1
736         massf1 = massf1 + atmass(i)
737      enddo
738
739   case('cma','cma interface')
740      ! copied from molbld.f without modification (except for error handling)
741      if (narg.ne.1) then
742         call env%error('not enough argument to constrain cma',source)
743         return
744      endif
745      if (key.eq.'cma interface') call cmaiface(nat,at,xyz)
746      nconstr = nconstr+1
747      atconstr(1,nconstr) = 0
748      if (trim(argv(1)).eq.'auto') then
749         massf1 = 0.0_wp
750         massf2 = 0.0_wp
751         do i = 1, nat
752            if (splitlist(i).eq.1) then
753               massf1 = massf1 + atmass(i)
754            else
755               massf2 = massf2 + atmass(i)
756            endif
757         enddo
758         call cmafrag(nat,at,xyz,ra,rb)
759         valconstr(nconstr) = rcma
760         write(env%unit,'("constraining fragment CMA to initial R(Ang.)=",f8.3)')&
761            valconstr(nconstr)*autoaa
762      else
763         if (getValue(env,trim(argv(1)),ddum)) then
764            valconstr(nconstr) = ddum*aatoau
765         else
766            call env%warning("Something went wrong in set_constr_ 'cma'.",source)
767            nconstr = nconstr-1
768            return
769         endif
770      endif
771
772!   case('x','y','z')
773!      idum = index('xyz',key)
774   case('z')
775      ! copied from molbld.f without modification (except for error handling)
776      if (narg.ne.1) then
777         call env%error('not enough argument to constrain z coordinate',source)
778         return
779      endif
780      zconstr = 1
781      nconstr = nconstr+1
782      atconstr(1,nconstr) = -1
783      if (getValue(env,trim(argv(1)),ddum)) then
784         valconstr(nconstr) = ddum*aatoau
785      else
786         call env%warning("Something went wrong in set_constr_ 'z'.",source)
787         nconstr = nconstr-1
788         return
789      endif
790      if (sum(abs(xyz(3,1:iatf1))).gt.1.d-3) then
791         call env%warning('z-coordinates of fragment 1 must be 0!',source)
792         nconstr = nconstr-1
793         return
794      endif
795      write(env%unit,'("constraining fragment 2 to Z=0 plane at (Ang.)=",f8.3)') &
796         valconstr(nconstr)*autoaa
797
798   end select
799end subroutine set_constr
800
801!! --------------------------------------------------------------[SAW1809]-
802!  this is the new version of the scan routine exploiting all features
803subroutine set_scan(env,key,val,nat,at,idMap,xyz)
804   use xtb_scanparam
805   use xtb_mctc_convert, only : aatoau
806   use xtb_mctc_constants, only : pi
807   implicit none
808   character(len=*), parameter :: source = 'userdata_scan'
809   type(TEnvironment), intent(inout) :: env
810   character(len=*),intent(in) :: key
811   character(len=*),intent(in) :: val
812   integer, intent(in) :: nat
813   integer, intent(in) :: at(nat)
814   type(TIdentityMap), intent(in) :: idMap
815   real(wp),intent(in) :: xyz(3,nat)
816
817   integer  :: idum
818   real(wp) :: ddum
819   logical  :: ldum
820
821   integer  :: i,ie
822   real(wp) :: start_value,end_value
823   integer  :: narg
824   character(len=p_str_length),dimension(p_arg_length) :: argv
825   character(len=:),allocatable :: temp
826
827   nscan = nscan+1 ! make a new scan
828
829   ie = index(val,semicolon)
830   if (ie.ne.0) then
831      temp = val(:ie-1)
832      idum = nconstr
833      call set_constr(env,key,temp,nat,at,idMap,xyz) ! generate a new constraint
834      scan_list(nscan)%iconstr = nconstr ! new generated
835      if (idum.eq.nconstr) then
836         call env%error('Failed to generate constraint',source)
837         return
838      end if
839      temp = val(1+ie:)
840   else
841      if (getValue(env,key,idum)) then
842         if (idum.gt.nconstr) then
843            call env%error('Constraint '''//key//''' is not defined',source)
844            return
845         endif
846         scan_list(nscan)%iconstr = idum
847      else
848         call env%error('Constraint '''//key//''' is invalid in this context',source)
849         return
850      endif
851      temp = val
852   endif
853
854   call parse(temp,comma,argv,narg)
855   if (verbose) then
856      do idum = 1, narg
857         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
858      enddo
859   endif
860
861   ! get back the constraint number which should be scanned
862   i = scan_list(nscan)%iconstr
863   ! instead of simply saving the kind of constraint while reading,
864   ! we used some really obvious integer code system depending on
865   ! the number of elements in the atconstr-Array.
866   ! Expected to break at some point, add a FIXME and wait for complaints...
867   idum = atconstr(1,i)
868   if (atconstr(2,i).gt.0) idum = 1
869   if (atconstr(3,i).gt.0) idum = 2
870   if (atconstr(4,i).gt.0) idum = 3
871
872   if (getValue(env,trim(argv(1)),ddum)) then
873      if (idum.le.1) then
874         start_value = ddum * aatoau
875      else
876         start_value = ddum * pi/180.0_wp
877      endif
878   else
879      call env%error('Invalid start value for scan',source)
880      return
881   endif
882   if (getValue(env,trim(argv(2)),ddum)) then
883      if (idum.le.1) then
884         end_value = ddum * aatoau
885      else
886         end_value = ddum * pi/180.0_wp
887      endif
888   else
889      call env%error('Invalid end value for scan',source)
890      return
891   endif
892   if (getValue(env,trim(argv(3)),idum)) then
893      scan_list(nscan)%nscan = idum
894      allocate( scan_list(nscan)%valscan(idum), source = 0.0_wp )
895   else
896      call env%error('Invalid step number for scan',source)
897      return
898   endif
899
900   do i = 1, scan_list(nscan)%nscan
901      scan_list(nscan)%valscan(i) = start_value + (end_value-start_value) &
902         &                       * real(i-1,wp)/real(scan_list(nscan)%nscan-1,wp)
903      if (verbose) write(env%unit,'(i5,1x,f12.8)') i,scan_list(nscan)%valscan(i)
904   enddo
905
906end subroutine set_scan
907
908subroutine set_wall(env,key,val,nat,at,idMap,xyz)
909   use xtb_mctc_convert, only : autoaa
910   use xtb_sphereparam
911   implicit none
912   character(len=*), parameter :: source = 'userdata_wall'
913   type(TEnvironment), intent(inout) :: env
914   character(len=*),intent(in) :: key
915   character(len=*),intent(in) :: val
916   integer, intent(in) :: nat
917   integer, intent(in) :: at(nat)
918   type(TIdentityMap), intent(in) :: idMap
919   real(wp),intent(in) :: xyz(3,nat)
920
921   integer  :: idum
922   real(wp) :: ddum,darray(3)
923   logical  :: ldum
924   integer  :: list(nat),nlist
925   integer  :: tlist(nat),ntlist
926   integer  :: iarg,iaxis
927   real(wp) :: radius,radii(3),center(3)
928
929   integer  :: narg
930   character(len=p_str_length),dimension(p_arg_length) :: argv
931
932   nlist = 0
933   ntlist = 0
934   list = 0
935   tlist = 0
936
937   radius = 0.0_wp
938   radii  = 0.0_wp
939   center = 0.0_wp
940
941   call parse(val,comma,argv,narg)
942   if (verbose) then
943      do idum = 1, narg
944         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
945      enddo
946   endif
947   select case(key)
948   case default ! ignore, don't even think about raising them
949   case('sphere')
950      if (narg.lt.2) then
951         call env%warning("Not enough arguments to set up a spherical wall",source)
952         return
953      endif
954   !  part 1: get the sphere radius
955      if (trim(argv(1)).eq.'auto') then
956         call get_sphere_radius(nat,at,xyz,center,radius,do_trafo=.false.)
957      else
958         if (getValue(env,trim(argv(1)),ddum)) then
959            radius = ddum
960            center = 0.0_wp
961         else
962            ! warning already generated by get_value
963            return ! something went wrong
964         endif
965      endif
966   !  part 2: get atoms
967      if (trim(argv(2)).eq.'all') then
968         call set_sphere_radius(radius,center)
969      else
970         do iarg = 2, narg
971            if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then
972               if (nlist+ntlist.gt.nat) then
973                  call env%warning("Too many atoms in list for spherical wall.",source)
974                  return ! something went wrong
975               endif
976               if (maxval(tlist(:ntlist)).gt.nat) then
977                  call env%warning("Attempted to wall in a non-existing atom",source)
978                  cycle ! skip crappy input
979               endif
980               list(nlist+1:nlist+ntlist) = tlist
981               nlist = nlist + ntlist
982            else
983               ! warning already generated by get_list_value
984               return ! something went wrong
985            endif
986         enddo
987         call set_sphere_radius(radius,center,nlist,list)
988      endif
989      write(env%unit,'("spherical wallpotenial with radius",'//&
990         '1x,f12.7,1x,"Å")') radius*autoaa
991
992   case('ellipsoid')
993      if (narg.lt.4) then
994         call env%warning("Not enough arguments to set up an ellipsoidal wall",source)
995         return
996      endif
997   !  part 1: get ellipsoid axis
998      if ((trim(argv(1)).eq.'auto').or.(trim(argv(2)).eq.'auto').or. &
999          (trim(argv(3)).eq.'auto')) then
1000         ! use get_ellipsoid_radius if you manage to implement it
1001         call get_sphere_radius(nat,at,xyz,center,radius,do_trafo=.false.)
1002      endif
1003      do iaxis = 1, 3
1004         if (trim(argv(iaxis)).eq.'auto') then
1005            radii(iaxis) = radius
1006         else
1007            if (getValue(env,trim(argv(iaxis)),ddum)) then
1008               radii(iaxis) = ddum
1009               center(iaxis) = 0.0_wp
1010            else
1011               ! warning already generated by get_value
1012               return ! something went wrong
1013            endif
1014         endif
1015      enddo
1016   !  part 2: get atoms
1017      if (trim(argv(4)).eq.'all') then
1018         call set_sphere_radius(radii,center)
1019      else
1020         do iarg = 4, narg
1021            if (getListValue(env,trim(argv(iarg)),tlist,ntlist)) then
1022               if (nlist+ntlist.gt.nat) then
1023                  call env%warning("Too many atoms in list for spherical wall.",source)
1024                  return ! something went wrong
1025               endif
1026               if (maxval(tlist(:ntlist)).gt.nat) then
1027                  call env%warning("Attempted to wall in a non-existing atom",source)
1028                  cycle ! skip crappy input
1029               endif
1030               list(nlist+1:nlist+ntlist) = tlist
1031               nlist = nlist + ntlist
1032            else
1033               ! warning already generated by get_list_value
1034               return ! something went wrong
1035            endif
1036         enddo
1037         call set_sphere_radius(radii,center,nlist,list)
1038      endif
1039      write(env%unit,'("ellipsoidal wallpotenial with radii",'//&
1040         '3(1x,f12.7,1x,"Å"))') radii*autoaa
1041
1042   end select
1043
1044end subroutine set_wall
1045
1046subroutine set_split(env,key,val,nat,at,idMap,xyz)
1047   use xtb_splitparam
1048   implicit none
1049   character(len=*), parameter :: source = 'userdata_split'
1050   type(TEnvironment), intent(inout) :: env
1051   character(len=*),intent(in) :: key
1052   character(len=*),intent(in) :: val
1053   integer, intent(in) :: nat
1054   integer, intent(in) :: at(nat)
1055   type(TIdentityMap), intent(in) :: idMap
1056   real(wp),intent(in) :: xyz(3,nat)
1057
1058   integer  :: idum
1059   real(wp) :: ddum
1060   logical  :: ldum
1061   integer  :: list(nat),nlist
1062   integer  :: i,j
1063
1064   integer  :: narg
1065   character(len=p_str_length),dimension(p_arg_length) :: argv
1066
1067   call parse(val,comma,argv,narg)
1068   if (verbose) then
1069      do idum = 1, narg
1070         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1071      enddo
1072   endif
1073   select case(key)
1074   case default ! ignore, don't even think about raising them
1075   case('fragment1')
1076      do i = 1, narg
1077         if (getListValue(env,trim(argv(i)),list,nlist)) then
1078            do j = 1, nlist
1079               splitlist(list(j)) = 1
1080            enddo
1081            iatf1 = iatf1 + nlist ! compatibility stuff, as usual
1082            iatf(1) = iatf(1) + nlist
1083         else
1084            return ! you screwed it, let's get out of here
1085         endif
1086      enddo
1087   case('fragment2')
1088      do i = 1, narg
1089         if (getListValue(env,trim(argv(i)),list,nlist)) then
1090            do j = 1, nlist
1091               splitlist(list(j)) = 2
1092            enddo
1093            iatf2 = iatf2 + nlist ! compatibility stuff, as usual
1094            iatf(2) = iatf(2) + nlist
1095         else
1096            return ! you screwed it, let's get out of here
1097         endif
1098      enddo
1099   case('fragment')
1100      if (getValue(env,trim(argv(1)),idum)) then
1101         if (idum.gt.nat) then
1102            call env%warning("rejecting fragment number greater than number of atoms",source)
1103            return ! doesn't really make sense, sorry, your problem
1104         endif
1105         do i = 2, narg
1106            if (getListValue(env,trim(argv(i)),list,nlist)) then
1107               do j = 1, nlist
1108                  splitlist(list(j)) = idum
1109               enddo
1110               iatf(idum) = iatf(idum) + nlist
1111               if (idum.eq.1) then ! legacy, don't remove or something breaks
1112                  iatf1 = iatf1 + nlist
1113               else if (idum.eq.2) then
1114                  iatf2 = iatf2 + nlist
1115               endif
1116            else
1117               return ! you screwed it, let's get out of here
1118            endif
1119         enddo
1120      endif
1121   end select
1122end subroutine set_split
1123
1124subroutine set_hess(env,key,val,nat,at,idMap,xyz)
1125   use xtb_splitparam
1126   implicit none
1127   character(len=*), parameter :: source = 'userdata_hess'
1128   type(TEnvironment), intent(inout) :: env
1129   character(len=*),intent(in) :: key
1130   character(len=*),intent(in) :: val
1131   integer, intent(in) :: nat
1132   integer, intent(in) :: at(nat)
1133   type(TIdentityMap), intent(in) :: idMap
1134   real(wp),intent(in) :: xyz(3,nat)
1135
1136   integer  :: idum
1137   real(wp) :: ddum
1138   logical  :: ldum
1139   integer  :: i,j
1140   integer, allocatable :: list(:)
1141
1142   integer  :: narg
1143   character(len=p_str_length),dimension(p_arg_length) :: argv
1144
1145   call parse(val,comma,argv,narg)
1146   if (verbose) then
1147      do idum = 1, narg
1148         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1149      enddo
1150   endif
1151   select case(key)
1152   case default ! ignore, don't even think about raising them
1153   case('element mass')
1154      if (mod(narg,2).ne.0) then
1155         call env%warning("Something went wrong in set_mass_ 'element'.",source)
1156      endif
1157      do i = 1, narg, 2
1158         j = i+1
1159         if (getValue(env,trim(argv(j)),ddum)) then
1160            if (idMap%has(argv(i))) then
1161               call idMap%set(atmass, argv(i), ddum)
1162               write(env%unit,'(a,a,a,1x,g0)') &
1163                  "mass of elements '",trim(argv(i)),"' changed to", ddum
1164            else if (getValue(env,trim(argv(i)),idum)) then
1165               where(at.eq.idum) atmass = ddum
1166               write(env%unit,'(a,i0,a,1x,g0)') &
1167                  "mass of elements with Z=",idum," changed to", ddum
1168            end if
1169         end if
1170      enddo
1171   case('modify mass','isotope')
1172      if (mod(narg,2).ne.0) then
1173         call env%warning("Something went wrong in set_mass_ 'modify'.",source)
1174      endif
1175      do i = 1, narg, 2
1176         j = i+1
1177         if (getValue(env,trim(argv(i)),idum).and.&
1178             getValue(env,trim(argv(j)),ddum)) then
1179            if (idum.gt.nat) then
1180               call env%warning('Attempted setting atom mass not present in system.',source)
1181               cycle
1182            endif
1183            atmass(idum) = ddum
1184            write(env%unit,'(a,1x,i0,1x,a,1x,g0)') &
1185               'mass of atom ',idum,' changed to',atmass(idum)
1186         endif
1187      enddo
1188   case('scale mass')
1189      if (mod(narg,2).ne.0) then
1190         call env%warning("Something went wrong in set_mass_ 'scale'.",source)
1191      endif
1192      do i = 1, narg, 2
1193         j = i+1
1194         if (getValue(env,trim(argv(i)),idum).and.&
1195             getValue(env,trim(argv(j)),ddum)) then
1196            if (idum.gt.nat) then
1197               call env%warning('Attempted scaling atom not present in system.',source)
1198               cycle
1199            endif
1200            atmass(idum) = atmass(idum)*ddum
1201            write(env%unit,'(a,1x,i0,1x,a,1x,g0)') &
1202               'mass of atom ',idum,' changed to',atmass(idum)
1203         endif
1204      enddo
1205   end select
1206
1207end subroutine set_hess
1208
1209subroutine set_reactor(env,key,val,nat,at,idMap,xyz)
1210   use xtb_type_atomlist
1211   use xtb_setparam
1212   implicit none
1213   character(len=*), parameter :: source = 'userdata_reactor'
1214   type(TEnvironment), intent(inout) :: env
1215   character(len=*),intent(in) :: key
1216   character(len=*),intent(in) :: val
1217   integer, intent(in) :: nat
1218   integer, intent(in) :: at(nat)
1219   type(TIdentityMap), intent(in) :: idMap
1220   real(wp),intent(in) :: xyz(3,nat)
1221
1222   type(TAtomList) :: atl
1223
1224   integer  :: idum
1225   real(wp) :: ddum
1226   logical  :: ldum
1227   integer  :: nlist
1228   integer, allocatable :: list(:)
1229   integer  :: i,j
1230
1231   integer  :: narg
1232   character(len=p_str_length),dimension(p_arg_length) :: argv
1233
1234   call parse(val,comma,argv,narg)
1235   if (verbose) then
1236      do idum = 1, narg
1237         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1238      enddo
1239   endif
1240   select case(key)
1241   case default ! ignore, don't even think about raising them
1242   case('atoms')
1243      call atl%new(val)
1244      if (atl%get_error()) then
1245         call env%warning('something is wrong in the reactor atom list',source)
1246         return
1247      endif
1248      if (reactset%nat > 0) call atl%add(reactset%atoms(:reactset%nat))
1249      call atl%to_list(list)
1250      reactset%atoms = list
1251      reactset%nat = size(list)
1252   end select
1253
1254end subroutine set_reactor
1255
1256subroutine set_path(env,key,val,nat,at,idMap,xyz)
1257   use xtb_type_atomlist
1258   use xtb_setparam
1259   implicit none
1260   character(len=*), parameter :: source = 'userdata_path'
1261   type(TEnvironment), intent(inout) :: env
1262   character(len=*),intent(in) :: key
1263   character(len=*),intent(in) :: val
1264   integer, intent(in) :: nat
1265   integer, intent(in) :: at(nat)
1266   type(TIdentityMap), intent(in) :: idMap
1267   real(wp),intent(in) :: xyz(3,nat)
1268
1269   type(TAtomList) :: atl
1270
1271   integer  :: idum
1272   real(wp) :: ddum
1273   logical  :: ldum
1274   integer  :: nlist
1275   integer, allocatable :: list(:)
1276   integer  :: i,j
1277
1278   integer  :: narg
1279   character(len=p_str_length),dimension(p_arg_length) :: argv
1280
1281   call parse(val,comma,argv,narg)
1282   if (verbose) then
1283      do idum = 1, narg
1284         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1285      enddo
1286   endif
1287   select case(key)
1288   case default ! ignore, don't even think about raising them
1289   case('atoms')
1290      call atl%new(val)
1291      if (atl%get_error()) then
1292         call env%warning('something is wrong in the bias path atom list',source)
1293         return
1294      endif
1295      if (pathset%nat > 0) call atl%add(pathset%atoms(:pathset%nat))
1296      call atl%to_list(list)
1297      pathset%atoms = list
1298      pathset%nat = size(list)
1299   end select
1300
1301end subroutine set_path
1302
1303subroutine set_metadyn(env,key,val,nat,at,idMap,xyz)
1304   use xtb_type_atomlist
1305   use xtb_fixparam
1306   implicit none
1307   character(len=*), parameter :: source = 'userdata_metadyn'
1308   type(TEnvironment), intent(inout) :: env
1309   character(len=*),intent(in) :: key
1310   character(len=*),intent(in) :: val
1311   integer, intent(in) :: nat
1312   integer, intent(in) :: at(nat)
1313   type(TIdentityMap), intent(in) :: idMap
1314   real(wp),intent(in) :: xyz(3,nat)
1315
1316   type(TAtomList) :: atl
1317
1318   integer  :: idum
1319   real(wp) :: ddum
1320   logical  :: ldum
1321   integer  :: nlist
1322   integer, allocatable :: list(:)
1323   integer  :: i,j,iat
1324
1325   integer  :: narg
1326   character(len=p_str_length),dimension(p_arg_length) :: argv
1327
1328   call parse(val,comma,argv,narg)
1329   if (verbose) then
1330      do idum = 1, narg
1331         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1332      enddo
1333   endif
1334   select case(key)
1335   case default ! ignore, don't even think about raising them
1336   case('bias atoms')
1337      call atl%new(val)
1338      if (atl%get_error()) then
1339         call env%warning('something is wrong in the metadynamic atom list',source)
1340         return
1341      endif
1342      if (rmsdset%nat > 0) call atl%add(rmsdset%atoms(:rmsdset%nat))
1343      call atl%to_list(list)
1344      rmsdset%atoms = list
1345      rmsdset%nat = size(list)
1346   case('bias elements')
1347      call atl%new
1348      do idum = 1, narg
1349         ! get element by symbol
1350         if (idMap%has(argv(idum))) then
1351            call idMap%get(list, argv(idum))
1352            if (allocated(list)) then
1353               call atl%add(list)
1354            else
1355               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
1356               cycle
1357            end if
1358         else
1359            ldum = getValue(env,trim(argv(idum)),iat)
1360            if (.not.ldum) cycle ! skip garbage input
1361            ! check for unreasonable input
1362            if (iat > 0) then
1363               call atl%add(at.eq.iat)
1364            else
1365               call env%warning("Unknown element: '"//trim(argv(idum))//"'",source)
1366               cycle
1367            endif
1368         endif
1369      enddo
1370      if (rmsdset%nat > 0) call atl%add(rmsdset%atoms(:rmsdset%nat))
1371      call atl%to_list(list)
1372      rmsdset%atoms = list
1373      rmsdset%nat = size(list)
1374   case('atoms')
1375      call atl%new(val)
1376      if (atl%get_error()) then
1377         call env%warning('something is wrong in the metadynamic atom list',source)
1378         return
1379      endif
1380      if (metaset%nat > 0) call atl%add(metaset%atoms(:metaset%nat))
1381      call atl%to_list(list)
1382      metaset%atoms = list
1383      metaset%nat = size(list)
1384   case('modify factor')
1385      if (mod(narg,2).ne.0) then
1386         call env%warning("Something went wrong in set_metadyn_ 'modify'.",source)
1387      endif
1388      do i = 1, narg, 2
1389         j = i+1
1390         if (getValue(env,trim(argv(i)),idum).and.&
1391             getValue(env,trim(argv(j)),ddum)) then
1392            if (idum.gt.metaset%maxsave) then
1393               call env%warning('Attempted using factor not present in system.',source)
1394               cycle
1395            endif
1396            metaset%factor(idum) = ddum
1397         endif
1398      enddo
1399   case('scale factor')
1400      if (mod(narg,2).ne.0) then
1401         call env%warning("Something went wrong in set_metadyn_ 'scale'.",source)
1402      endif
1403      do i = 1, narg, 2
1404         j = i+1
1405         if (getValue(env,trim(argv(i)),idum).and.&
1406             getValue(env,trim(argv(j)),ddum)) then
1407            if (idum.gt.metaset%maxsave) then
1408               call env%warning('Attempted scaling factor not present in system.',source)
1409               cycle
1410            endif
1411            metaset%factor(idum) = metaset%factor(idum)*ddum
1412         endif
1413      enddo
1414   end select
1415
1416end subroutine set_metadyn
1417
1418subroutine set_freeze(env,key,val,nat,at,idMap,xyz)
1419   implicit none
1420   character(len=*), parameter :: source = 'userdata_freeze'
1421   type(TEnvironment), intent(inout) :: env
1422   character(len=*),intent(in) :: key
1423   character(len=*),intent(in) :: val
1424   integer, intent(in) :: nat
1425   integer, intent(in) :: at(nat)
1426   type(TIdentityMap), intent(in) :: idMap
1427   real(wp),intent(in) :: xyz(3,nat)
1428
1429   integer  :: idum
1430   real(wp) :: ddum
1431   logical  :: ldum
1432   integer  :: list(nat),nlist
1433   integer  :: i,j
1434
1435   integer  :: narg
1436   character(len=p_str_length),dimension(p_arg_length) :: argv
1437
1438   call parse(val,comma,argv,narg)
1439   if (verbose) then
1440      do idum = 1, narg
1441         write(env%unit,'("-->",1x,i0,":",1x,a)') idum, trim(argv(idum))
1442      enddo
1443   endif
1444   select case(key)
1445   case default ! ignore, don't even think about raising them
1446   case('hessf')
1447   case('hessa')
1448   end select
1449
1450end subroutine set_freeze
1451
1452subroutine set_legacy(env,key,val,nat,at,idMap,xyz)
1453   implicit none
1454   character(len=*), parameter :: source = 'userdata_legacy'
1455   type(TEnvironment), intent(inout) :: env
1456   character(len=*),intent(in) :: key
1457   character(len=*),intent(in) :: val
1458   integer, intent(in) :: nat
1459   integer, intent(in) :: at(nat)
1460   type(TIdentityMap), intent(in) :: idMap
1461   real(wp),intent(in) :: xyz(3,nat)
1462   integer  :: err
1463   integer  :: idum
1464   real(wp) :: ddum
1465   logical  :: ldum
1466   select case(key)
1467   case default ! complaining about unknown keywords should already be done
1468      continue  ! so we do nothing here
1469   case('hessf');       call set_fix(env,'freeze',val,nat,at,idMap,xyz)
1470!   case('hessa');       call set_frozh('hessa',val)
1471   case('fragment1'); call set_split(env,'fragment1',val,nat,at,idMap,xyz)
1472   case('fragment2'); call set_split(env,'fragment1',val,nat,at,idMap,xyz)
1473   case('constrxyz'); call set_fix(env,'atoms',val,nat,at,idMap,xyz)
1474!   case('constrainel')
1475!   case('constrain')
1476!   case('scan')
1477   case('ellips'); call set_wall(env,'ellipsoid',val,nat,at,idMap,xyz)
1478   case('sphere'); call set_wall(env,'sphere',val,nat,at,idMap,xyz)
1479   case('fix'); call set_fix(env,'atoms',val,nat,at,idMap,xyz)
1480   case('atomlist+'); call set_metadyn(env,'atoms',val,nat,at,idMap,xyz)
1481   end select
1482
1483end subroutine set_legacy
1484
1485end module xtb_constrain_param
1486