CP2K 8.1 compact CDFT weight-evaluator source excerpts
runtime revision: git:0b61f2f
container source HEAD: 0bfbd41a7402152da7b86a103d5d95a4184367d9
One-based line numbers refer to the exact embedded /cp2k/src files.

===== qs_cdft_utils.F:100-110 =====
   100     SUBROUTINE becke_constraint_init(qs_env)
   101        TYPE(qs_environment_type), POINTER                 :: qs_env
   102  
   103        CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_init'
   104  
   105        CHARACTER(len=2)                                   :: element_symbol
   106        INTEGER :: atom_a, bounds(2), handle, i, iatom, iex, igroup, ikind, ip, ithread, iw, j, &
   107           jatom, katom, natom, nkind, npme, nthread, numexp, unit_nr
   108        INTEGER, DIMENSION(2, 3)                           :: bo
   109        INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores, stride
   110        LOGICAL                                            :: build, in_memory, mpi_io

===== qs_cdft_utils.F:167-275 =====
   167        IF (becke_control%adjust) THEN
   168           IF (.NOT. ASSOCIATED(becke_control%radii)) THEN
   169              CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
   170              IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%radii_tmp)) &
   171                 CALL cp_abort(__LOCATION__, &
   172                               "Length of keyword BECKE_CONSTRAINT\ATOMIC_RADII does not "// &
   173                               "match number of atomic kinds in the input coordinate file.")
   174              ALLOCATE (becke_control%radii(SIZE(atomic_kind_set)))
   175              becke_control%radii(:) = becke_control%radii_tmp(:)
   176              DEALLOCATE (becke_control%radii_tmp)
   177           END IF
   178        END IF
   179        ! Setup cutoff scheme
   180        IF (.NOT. ASSOCIATED(becke_control%cutoffs)) THEN
   181           CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
   182           ALLOCATE (becke_control%cutoffs(natom))
   183           SELECT CASE (becke_control%cutoff_type)
   184           CASE (becke_cutoff_global)
   185              becke_control%cutoffs(:) = becke_control%rglobal
   186           CASE (becke_cutoff_element)
   187              IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%cutoffs_tmp)) &
   188                 CALL cp_abort(__LOCATION__, &
   189                               "Length of keyword BECKE_CONSTRAINT\ELEMENT_CUTOFFS does not "// &
   190                               "match number of atomic kinds in the input coordinate file.")
   191              DO ikind = 1, SIZE(atomic_kind_set)
   192                 CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
   193                 DO iatom = 1, katom
   194                    atom_a = atom_list(iatom)
   195                    becke_control%cutoffs(atom_a) = becke_control%cutoffs_tmp(ikind)
   196                 END DO
   197              END DO
   198              DEALLOCATE (becke_control%cutoffs_tmp)
   199           END SELECT
   200        END IF
   201        ! Zero weight functions
   202        DO igroup = 1, SIZE(group)
   203           group(igroup)%weight%pw%cr3d = 0.0_dp
   204        END DO
   205        IF (cdft_control%atomic_charges) THEN
   206           DO iatom = 1, cdft_control%natoms
   207              cdft_control%charge(iatom)%pw%cr3d = 0.0_dp
   208           END DO
   209        END IF
   210        ! Allocate storage for cell adjustment coefficients and needed distance vectors
   211        build = .FALSE.
   212        IF (becke_control%adjust .AND. .NOT. ASSOCIATED(becke_control%aij)) THEN
   213           ALLOCATE (becke_control%aij(natom, natom))
   214           build = .TRUE.
   215        END IF
   216        IF (becke_control%vector_buffer%store_vectors) THEN
   217           ALLOCATE (becke_control%vector_buffer%distances(natom))
   218           ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
   219           IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
   220           ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
   221        END IF
   222        ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
   223        ! Calculate pairwise distances between each atom pair
   224        DO i = 1, 3
   225           cell_v(i) = cell%hmat(i, i)
   226        END DO
   227        DO iatom = 1, natom - 1
   228           DO jatom = iatom + 1, natom
   229              r = particle_set(iatom)%r
   230              r1 = particle_set(jatom)%r
   231              DO i = 1, 3
   232                 r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
   233                 r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
   234              END DO
   235              dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
   236              ! Store pbc corrected position and pairwise distance vectors for later reuse
   237              IF (becke_control%vector_buffer%store_vectors) THEN
   238                 becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
   239                 IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
   240                 IF (in_memory) THEN
   241                    becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
   242                    becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
   243                 END IF
   244              END IF
   245              becke_control%vector_buffer%R12(iatom, jatom) = SQRT(DOT_PRODUCT(dist_vec, dist_vec))
   246              becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
   247              ! Set up heteronuclear cell partitioning using user defined radii
   248              IF (build) THEN
   249                 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=ikind)
   250                 ircov = becke_control%radii(ikind)
   251                 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, kind_number=ikind)
   252                 jrcov = becke_control%radii(ikind)
   253                 IF (ircov .NE. jrcov) THEN
   254                    chi = ircov/jrcov
   255                    uij = (chi - 1.0_dp)/(chi + 1.0_dp)
   256                    becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp)
   257                    IF (becke_control%aij(iatom, jatom) .GT. 0.5_dp) THEN
   258                       becke_control%aij(iatom, jatom) = 0.5_dp
   259                    ELSE IF (becke_control%aij(iatom, jatom) .LT. -0.5_dp) THEN
   260                       becke_control%aij(iatom, jatom) = -0.5_dp
   261                    END IF
   262                 ELSE
   263                    becke_control%aij(iatom, jatom) = 0.0_dp
   264                 END IF
   265                 ! Note change of sign
   266                 becke_control%aij(jatom, iatom) = -becke_control%aij(iatom, jatom)
   267              END IF
   268           END DO
   269        END DO
   270        ! Dump some additional information about the calculation
   271        IF (cdft_control%first_iteration) THEN
   272           IF (iw > 0) THEN
   273              WRITE (iw, '(/,T3,A)') &
   274                 '----------------------- Becke atomic parameters ------------------------'
   275              IF (becke_control%adjust) THEN

===== qs_cdft_utils.F:534-610 =====
   534     SUBROUTINE read_constraint_definitions(cdft_control, cdft_control_section)
   535  
   536        TYPE(cdft_control_type), INTENT(INOUT)             :: cdft_control
   537        TYPE(section_vals_type), INTENT(INOUT), POINTER    :: cdft_control_section
   538  
   539        INTEGER                                            :: i, j, jj, k, n_rep, natoms, nvar, &
   540                                                              tot_natoms
   541        INTEGER, DIMENSION(:), POINTER                     :: atomlist, dummylist, tmplist
   542        LOGICAL                                            :: exists, is_duplicate
   543        REAL(KIND=dp), DIMENSION(:), POINTER               :: rtmplist
   544        TYPE(section_vals_type), POINTER                   :: group_section
   545  
   546        NULLIFY (tmplist, rtmplist, atomlist, dummylist, group_section)
   547  
   548        group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP")
   549        CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists)
   550        IF (.NOT. exists) CPABORT("Section ATOM_GROUP is missing.")
   551        ALLOCATE (cdft_control%group(nvar))
   552        tot_natoms = 0
   553        ! Parse all ATOM_GROUP sections
   554        DO k = 1, nvar
   555           ! First determine how much storage is needed
   556           natoms = 0
   557           CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, n_rep_val=n_rep)
   558           DO j = 1, n_rep
   559              CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist)
   560              IF (SIZE(tmplist) < 1) &
   561                 CPABORT("Each ATOM_GROUP must contain at least 1 atom.")
   562              natoms = natoms + SIZE(tmplist)
   563           END DO
   564           ALLOCATE (cdft_control%group(k)%atoms(natoms))
   565           ALLOCATE (cdft_control%group(k)%coeff(natoms))
   566           NULLIFY (cdft_control%group(k)%weight%pw)
   567           NULLIFY (cdft_control%group(k)%gradients)
   568           NULLIFY (cdft_control%group(k)%integrated)
   569           tot_natoms = tot_natoms + natoms
   570           ! Now parse
   571           jj = 0
   572           DO j = 1, n_rep
   573              CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist)
   574              DO i = 1, SIZE(tmplist)
   575                 jj = jj + 1
   576                 cdft_control%group(k)%atoms(jj) = tmplist(i)
   577              END DO
   578           END DO
   579           CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, n_rep_val=n_rep)
   580           jj = 0
   581           DO j = 1, n_rep
   582              CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, i_rep_val=j, r_vals=rtmplist)
   583              DO i = 1, SIZE(rtmplist)
   584                 jj = jj + 1
   585                 IF (jj > natoms) CPABORT("Length of keywords ATOMS and COEFF must match.")
   586                 IF (ABS(rtmplist(i)) /= 1.0_dp) CPABORT("Keyword COEFF accepts only values +/-1.0")
   587                 cdft_control%group(k)%coeff(jj) = rtmplist(i)
   588              END DO
   589           END DO
   590           IF (jj < natoms) CPABORT("Length of keywords ATOMS and COEFF must match.")
   591           CALL section_vals_val_get(group_section, "CONSTRAINT_TYPE", i_rep_section=k, &
   592                                     i_val=cdft_control%group(k)%constraint_type)
   593           CALL section_vals_val_get(group_section, "FRAGMENT_CONSTRAINT", i_rep_section=k, &
   594                                     l_val=cdft_control%group(k)%is_fragment_constraint)
   595           IF (cdft_control%group(k)%is_fragment_constraint) cdft_control%fragment_density = .TRUE.
   596        END DO
   597        ! Create a list containing all constraint atoms
   598        ALLOCATE (atomlist(tot_natoms))
   599        atomlist = -1
   600        jj = 0
   601        DO k = 1, nvar
   602           DO j = 1, SIZE(cdft_control%group(k)%atoms)
   603              is_duplicate = .FALSE.
   604              DO i = 1, jj + 1
   605                 IF (cdft_control%group(k)%atoms(j) == atomlist(i)) THEN
   606                    is_duplicate = .TRUE.
   607                    EXIT
   608                 END IF
   609              END DO
   610              IF (.NOT. is_duplicate) THEN

===== qs_cdft_utils.F:690-717 =====
   690        END DO
   691        CALL section_vals_val_get(cdft_control_section, "TARGET", r_vals=rtmplist)
   692        IF (SIZE(rtmplist) /= nvar) &
   693           CALL cp_abort(__LOCATION__, &
   694                         "The length of keyword TARGET is incorrect. "// &
   695                         "Expected "//TRIM(ADJUSTL(cp_to_string(nvar)))// &
   696                         " value(s), got "// &
   697                         TRIM(ADJUSTL(cp_to_string(SIZE(rtmplist))))//" value(s).")
   698        DO j = 1, nvar
   699           cdft_control%target(j) = rtmplist(j)
   700        END DO
   701        ! Read fragment constraint definitions
   702        IF (cdft_control%fragment_density) THEN
   703           CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_FILE_NAME", &
   704                                     c_val=cdft_control%fragment_a_fname)
   705           CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_FILE_NAME", &
   706                                     c_val=cdft_control%fragment_b_fname)
   707           CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_SPIN_FILE", &
   708                                     c_val=cdft_control%fragment_a_spin_fname)
   709           CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_SPIN_FILE", &
   710                                     c_val=cdft_control%fragment_b_spin_fname)
   711           CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_A", &
   712                                     l_val=cdft_control%flip_fragment(1))
   713           CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_B", &
   714                                     l_val=cdft_control%flip_fragment(2))
   715        END IF
   716  
   717     END SUBROUTINE read_constraint_definitions

===== qs_cdft_utils.F:966-986 =====
   966     SUBROUTINE hirshfeld_constraint_init(qs_env)
   967        TYPE(qs_environment_type), POINTER                 :: qs_env
   968  
   969        CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_init'
   970  
   971        CHARACTER(len=2)                                   :: element_symbol
   972        INTEGER                                            :: handle, iat, iatom, igroup, ikind, ip, &
   973                                                              iw, natom, nkind
   974        INTEGER, DIMENSION(:), POINTER                     :: atom_list
   975        REAL(KIND=dp)                                      :: zeff
   976        REAL(KIND=dp), DIMENSION(:), POINTER               :: radii_list
   977        TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
   978        TYPE(atomic_kind_type), POINTER                    :: atomic_kind
   979        TYPE(cdft_control_type), POINTER                   :: cdft_control
   980        TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
   981        TYPE(cp_logger_type), POINTER                      :: logger
   982        TYPE(dft_control_type), POINTER                    :: dft_control
   983        TYPE(hirshfeld_constraint_type), POINTER           :: hirshfeld_control
   984        TYPE(hirshfeld_type), POINTER                      :: hirshfeld_env
   985        TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
   986        TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set

===== qs_cdft_utils.F:1020-1086 =====
  1020                    radii_list(ikind) = cp_unit_from_cp2k(hirshfeld_control%radii(ikind), "angstrom")
  1021                 END IF
  1022              END DO
  1023           END IF
  1024           ! radius/radii_list parameters are optional for shape_function_density
  1025           CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, &
  1026                                      radius=hirshfeld_control%radius, &
  1027                                      radii_list=radii_list)
  1028           IF (ASSOCIATED(radii_list)) DEALLOCATE (radii_list)
  1029        END IF
  1030  
  1031        ! Atomic reference charges (Mulliken not supported)
  1032        IF (.NOT. ASSOCIATED(hirshfeld_env%charges)) THEN
  1033           CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
  1034                           nkind=nkind, natom=natom)
  1035           ALLOCATE (hirshfeld_env%charges(natom))
  1036           DO ikind = 1, nkind
  1037              CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
  1038              atomic_kind => atomic_kind_set(ikind)
  1039              CALL get_atomic_kind(atomic_kind, atom_list=atom_list)
  1040              DO iat = 1, SIZE(atom_list)
  1041                 iatom = atom_list(iat)
  1042                 hirshfeld_env%charges(iatom) = zeff
  1043              END DO
  1044           END DO
  1045        END IF
  1046  
  1047        ! Print some additional information about the calculation on the first iteration
  1048        IF (cdft_control%first_iteration) THEN
  1049           IF (iw > 0) THEN
  1050              group => cdft_control%group
  1051              CALL get_qs_env(qs_env, particle_set=particle_set)
  1052              IF (ASSOCIATED(hirshfeld_control%radii)) THEN
  1053                 WRITE (iw, '(T3,A)') &
  1054                    'Atom   Element  Gaussian radius (angstrom)'
  1055                 DO iatom = 1, natom
  1056                    CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
  1057                    WRITE (iw, "(i7,T15,A2,T37,F8.3)") &
  1058                       iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(hirshfeld_control%radii(iatom), "angstrom")
  1059                 END DO
  1060                 WRITE (iw, '(T3,A)') &
  1061                    '------------------------------------------------------------------------'
  1062              END IF
  1063              WRITE (iw, '(/,T3,A,T60)') &
  1064                 '----------------------- CDFT group definitions -------------------------'
  1065              DO igroup = 1, SIZE(group)
  1066                 IF (igroup > 1) WRITE (iw, '(T3,A)') ' '
  1067                 WRITE (iw, '(T5,A,I5,A,I5)') &
  1068                    'Atomic group', igroup, ' of ', SIZE(group)
  1069                 WRITE (iw, '(T5,A)') 'Atom  Element  Coefficient'
  1070                 DO ip = 1, SIZE(group(igroup)%atoms)
  1071                    iatom = group(igroup)%atoms(ip)
  1072                    CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
  1073                    WRITE (iw, '(i8,T16,A2,T23,F8.3)') iatom, ADJUSTR(element_symbol), group(igroup)%coeff(ip)
  1074                 END DO
  1075              END DO
  1076              WRITE (iw, '(T3,A)') &
  1077                 '------------------------------------------------------------------------'
  1078           END IF
  1079           cdft_control%first_iteration = .FALSE.
  1080        END IF
  1081  
  1082        ! Radii no longer needed
  1083        IF (ASSOCIATED(hirshfeld_control%radii)) DEALLOCATE (hirshfeld_control%radii)
  1084        CALL timestop(handle)
  1085  
  1086     END SUBROUTINE hirshfeld_constraint_init

===== qs_cdft_utils.F:1095-1130 =====
  1095     SUBROUTINE cdft_constraint_print(qs_env, electronic_charge)
  1096        TYPE(qs_environment_type), POINTER                 :: qs_env
  1097        REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: electronic_charge
  1098  
  1099        CHARACTER(len=2)                                   :: element_symbol
  1100        INTEGER                                            :: iatom, ikind, iw, jatom
  1101        REAL(kind=dp)                                      :: tc(2), zeff
  1102        TYPE(cdft_control_type), POINTER                   :: cdft_control
  1103        TYPE(cp_logger_type), POINTER                      :: logger
  1104        TYPE(dft_control_type), POINTER                    :: dft_control
  1105        TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
  1106        TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
  1107        TYPE(section_vals_type), POINTER                   :: cdft_constraint_section
  1108  
  1109        NULLIFY (cdft_constraint_section, logger, particle_set, dft_control, qs_kind_set)
  1110        logger => cp_get_default_logger()
  1111  
  1112        CALL get_qs_env(qs_env, &
  1113                        particle_set=particle_set, &
  1114                        dft_control=dft_control, &
  1115                        qs_kind_set=qs_kind_set)
  1116        CPASSERT(ASSOCIATED(qs_kind_set))
  1117  
  1118        cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
  1119        iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
  1120        cdft_control => dft_control%qs_control%cdft_control
  1121  
  1122        ! Print constraint information
  1123        CALL qs_scf_cdft_constraint_info(iw, cdft_control)
  1124  
  1125        ! Print weight function(s) to cube file(s) whenever weight is (re)built
  1126        IF (cdft_control%print_weight .AND. cdft_control%need_pot) &
  1127           CALL cdft_print_weight_function(qs_env)
  1128  
  1129        ! Print atomic CDFT charges
  1130        IF (iw > 0 .AND. cdft_control%atomic_charges) THEN

===== qs_cdft_utils.F:1257-1301 =====
  1257     SUBROUTINE cdft_print_weight_function(qs_env)
  1258        TYPE(qs_environment_type), POINTER                 :: qs_env
  1259  
  1260        CHARACTER(LEN=default_path_length)                 :: middle_name
  1261        INTEGER                                            :: igroup, unit_nr
  1262        LOGICAL                                            :: mpi_io
  1263        TYPE(cdft_control_type), POINTER                   :: cdft_control
  1264        TYPE(cp_logger_type), POINTER                      :: logger
  1265        TYPE(cp_para_env_type), POINTER                    :: para_env
  1266        TYPE(dft_control_type), POINTER                    :: dft_control
  1267        TYPE(particle_list_type), POINTER                  :: particles
  1268        TYPE(qs_subsys_type), POINTER                      :: subsys
  1269        TYPE(section_vals_type), POINTER                   :: cdft_constraint_section
  1270  
  1271        NULLIFY (cdft_constraint_section, logger, particles, dft_control, &
  1272                 para_env, subsys, cdft_control)
  1273        logger => cp_get_default_logger()
  1274  
  1275        CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control)
  1276        CALL qs_subsys_get(subsys, particles=particles)
  1277        cdft_control => dft_control%qs_control%cdft_control
  1278        cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
  1279  
  1280        DO igroup = 1, SIZE(cdft_control%group)
  1281           mpi_io = .TRUE.
  1282           middle_name = "cdft_weight_"//TRIM(ADJUSTL(cp_to_string(igroup)))
  1283           unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", &
  1284                                          middle_name=middle_name, &
  1285                                          extension=".cube", file_position="REWIND", &
  1286                                          log_filename=.FALSE., mpi_io=mpi_io)
  1287           ! Note PROGRAM_RUN_INFO section neeeds to be active!
  1288           IF (para_env%ionode .AND. unit_nr .LT. 1) &
  1289              CALL cp_abort(__LOCATION__, &
  1290                            "Please turn on PROGRAM_RUN_INFO to print CDFT weight function.")
  1291  
  1292           CALL cp_pw_to_cube(cdft_control%group(igroup)%weight%pw, &
  1293                              unit_nr, &
  1294                              "CDFT Weight Function", &
  1295                              particles=particles, &
  1296                              stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"), &
  1297                              mpi_io=mpi_io)
  1298           CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
  1299        END DO
  1300  
  1301     END SUBROUTINE cdft_print_weight_function

===== qs_cdft_methods.F:142-160 =====
   142     SUBROUTINE becke_constraint_low(qs_env, just_gradients)
   143        TYPE(qs_environment_type), POINTER                 :: qs_env
   144        LOGICAL, OPTIONAL                                  :: just_gradients
   145  
   146        CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_low'
   147  
   148        INTEGER                                            :: handle, i, iatom, igroup, ind(3), ip, j, &
   149                                                              jatom, jp, k, natom, np(3), nskipped
   150        INTEGER, ALLOCATABLE, DIMENSION(:)                 :: catom
   151        INTEGER, DIMENSION(2, 3)                           :: bo, bo_conf
   152        LOGICAL                                            :: in_memory, my_just_gradients
   153        LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: is_constraint, skip_me
   154        LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: atom_in_group
   155        REAL(kind=dp)                                      :: dist1, dist2, dmyexp, dvol, eps_cavity, &
   156                                                              my1, my1_homo, myexp, sum_cell_f_all, &
   157                                                              th, tmp_const
   158        REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: cell_functions, ds_dR_i, ds_dR_j, &
   159                                                              sum_cell_f_group
   160        REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: d_sum_Pm_dR, dP_i_dRi

===== qs_cdft_methods.F:293-320 =====
   293        ! Build constraint
   294        DO k = bo(1, 1), bo(2, 1)
   295           DO j = bo(1, 2), bo(2, 2)
   296              DO i = bo(1, 3), bo(2, 3)
   297                 ! If the grid point is too far from all constraint atoms and cavity confinement is active,
   298                 ! we can skip this grid point as it does not contribute to the weight or gradients
   299                 IF (becke_control%cavity_confine) THEN
   300                    IF (becke_control%cavity%pw%cr3d(k, j, i) < eps_cavity) CYCLE
   301                 END IF
   302                 ind = (/k, j, i/)
   303                 grid_p(1) = k*dr(1) + shift(1)
   304                 grid_p(2) = j*dr(2) + shift(2)
   305                 grid_p(3) = i*dr(3) + shift(3)
   306                 nskipped = 0
   307                 cell_functions = 1.0_dp
   308                 skip_me = .FALSE.
   309                 IF (becke_control%vector_buffer%store_vectors) becke_control%vector_buffer%distances = 0.0_dp
   310                 IF (in_memory) THEN
   311                    d_sum_Pm_dR = 0.0_dp
   312                    DO igroup = 1, SIZE(group)
   313                       group(igroup)%d_sum_const_dR = 0.0_dp
   314                    END DO
   315                    dP_i_dRi = 0.0_dp
   316                 END IF
   317                 ! Iterate over all atoms in the system
   318                 DO iatom = 1, natom
   319                    IF (skip_me(iatom)) THEN
   320                       cell_functions(iatom) = 0.0_dp

===== qs_cdft_methods.F:500-540 =====
   500                       skip_me(iatom) = .TRUE.
   501                       IF (becke_control%should_skip) THEN
   502                          IF (is_constraint(iatom)) nskipped = nskipped + 1
   503                          IF (nskipped == cdft_control%natoms) THEN
   504                             IF (in_memory) THEN
   505                                IF (becke_control%cavity_confine) THEN
   506                                   becke_control%cavity%pw%cr3d(k, j, i) = 0.0_dp
   507                                END IF
   508                             END IF
   509                             EXIT
   510                          END IF
   511                       END IF
   512                    END IF
   513                 END DO !iatom
   514                 IF (nskipped == cdft_control%natoms) CYCLE
   515                 ! Sum up cell functions
   516                 sum_cell_f_group = 0.0_dp
   517                 DO igroup = 1, SIZE(group)
   518                    DO ip = 1, SIZE(group(igroup)%atoms)
   519                       sum_cell_f_group(igroup) = sum_cell_f_group(igroup) + group(igroup)%coeff(ip)* &
   520                                                  cell_functions(group(igroup)%atoms(ip))
   521                    END DO
   522                 END DO
   523                 sum_cell_f_all = 0.0_dp
   524                 DO ip = 1, natom
   525                    sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
   526                 END DO
   527                 ! Gradients at (k,j,i)
   528                 IF (in_memory .AND. ABS(sum_cell_f_all) .GT. 0.0_dp) THEN
   529                    DO igroup = 1, SIZE(group)
   530                       DO iatom = 1, natom
   531                          group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
   532                             group(igroup)%d_sum_const_dR(1:3, iatom)/sum_cell_f_all - sum_cell_f_group(igroup)* &
   533                             d_sum_Pm_dR(1:3, iatom)/(sum_cell_f_all**2)
   534                       END DO
   535                    END DO
   536                 END IF
   537                 ! Weight function(s) at (k,j,i)
   538                 IF (.NOT. my_just_gradients .AND. ABS(sum_cell_f_all) .GT. 0.000001) THEN
   539                    DO igroup = 1, SIZE(group)
   540                       group(igroup)%weight%pw%cr3d(k, j, i) = sum_cell_f_group(igroup)/sum_cell_f_all

===== qs_cdft_methods.F:772-795 =====
   772     SUBROUTINE cdft_constraint_integrate(qs_env)
   773        TYPE(qs_environment_type), POINTER                 :: qs_env
   774  
   775        CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_integrate'
   776  
   777        INTEGER                                            :: handle, i, iatom, igroup, ikind, ivar, &
   778                                                              iw, jatom, natom, nvar
   779        LOGICAL                                            :: is_becke, paw_atom
   780        REAL(kind=dp)                                      :: dvol, eps_cavity, sign
   781        REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: dE, strength, target_val
   782        REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: electronic_charge, gapw_offset
   783        TYPE(becke_constraint_type), POINTER               :: becke_control
   784        TYPE(cdft_control_type), POINTER                   :: cdft_control
   785        TYPE(cdft_group_type), DIMENSION(:), POINTER       :: group
   786        TYPE(cp_logger_type), POINTER                      :: logger
   787        TYPE(cp_para_env_type), POINTER                    :: para_env
   788        TYPE(dft_control_type), POINTER                    :: dft_control
   789        TYPE(mpole_rho_atom), DIMENSION(:), POINTER        :: mp_rho
   790        TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
   791        TYPE(pw_p_type), DIMENSION(:), POINTER             :: charge, rho_r
   792        TYPE(qs_energy_type), POINTER                      :: energy
   793        TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
   794        TYPE(qs_rho_type), POINTER                         :: rho
   795        TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole

===== qs_cdft_methods.F:817-865 =====
   817        IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
   818           CPABORT("Becke control has not been allocated.")
   819        group => cdft_control%group
   820        ! Initialize
   821        nvar = SIZE(cdft_control%target)
   822        ALLOCATE (strength(nvar))
   823        ALLOCATE (target_val(nvar))
   824        ALLOCATE (dE(nvar))
   825        strength(:) = cdft_control%strength(:)
   826        target_val(:) = cdft_control%target(:)
   827        dE = 0.0_dp
   828        dvol = group(1)%weight%pw%pw_grid%dvol
   829        IF (cdft_control%atomic_charges) THEN
   830           charge => cdft_control%charge
   831           ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
   832           electronic_charge = 0.0_dp
   833        END IF
   834        ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
   835        DO i = 1, dft_control%nspins
   836           DO igroup = 1, SIZE(group)
   837              SELECT CASE (group(igroup)%constraint_type)
   838              CASE (cdft_charge_constraint)
   839                 sign = 1.0_dp
   840              CASE (cdft_magnetization_constraint)
   841                 IF (i == 1) THEN
   842                    sign = 1.0_dp
   843                 ELSE
   844                    sign = -1.0_dp
   845                 END IF
   846              CASE (cdft_alpha_constraint)
   847                 sign = 1.0_dp
   848                 IF (i == 2) CYCLE
   849              CASE (cdft_beta_constraint)
   850                 sign = 1.0_dp
   851                 IF (i == 1) CYCLE
   852              CASE DEFAULT
   853                 CPABORT("Unknown constraint type.")
   854              END SELECT
   855              IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
   856                 ! With external control, we can use cavity_mat as a mask to kahan sum
   857                 eps_cavity = becke_control%eps_cavity
   858                 IF (igroup /= 1) &
   859                    CALL cp_abort(__LOCATION__, &
   860                                  "Multiple constraints not yet supported by parallel mixed calculations.")
   861                 dE(igroup) = dE(igroup) + sign*accurate_dot_product(group(igroup)%weight%pw%cr3d, rho_r(i)%pw%cr3d, &
   862                                                                     becke_control%cavity_mat, eps_cavity)*dvol
   863              ELSE
   864                 dE(igroup) = dE(igroup) + sign*accurate_sum(group(igroup)%weight%pw%cr3d*rho_r(i)%pw%cr3d)*dvol
   865              END IF

===== qs_cdft_methods.F:935-954 =====
   935                 dE(ivar) = dE(ivar) + gapw_offset(ivar, i)
   936              END DO
   937           END DO
   938           DEALLOCATE (gapw_offset)
   939        END IF
   940        ! Update constraint value and energy
   941        cdft_control%value(:) = dE(:)
   942        energy%cdft = 0.0_dp
   943        DO ivar = 1, nvar
   944           energy%cdft = energy%cdft + (dE(ivar) - target_val(ivar))*strength(ivar)
   945        END DO
   946        ! Print constraint info and atomic CDFT charges
   947        CALL cdft_constraint_print(qs_env, electronic_charge)
   948        ! Deallocate tmp storage
   949        DEALLOCATE (dE, strength, target_val)
   950        IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
   951        CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
   952        CALL timestop(handle)
   953  
   954     END SUBROUTINE cdft_constraint_integrate

===== qs_cdft_methods.F:1153-1170 =====
  1153     SUBROUTINE hirshfeld_constraint_low(qs_env)
  1154        TYPE(qs_environment_type), POINTER                 :: qs_env
  1155  
  1156        CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_low'
  1157  
  1158        INTEGER                                            :: atom_a, handle, i, iatom, iex, igroup, &
  1159                                                              ikind, ithread, j, natom, npme, &
  1160                                                              nthread, numexp, subpatch_pattern
  1161        INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
  1162        LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: compute_charge, is_constraint
  1163        REAL(kind=dp)                                      :: alpha, coef, dvol, eps_rho_rspace, &
  1164                                                              radius, radius_constr
  1165        REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: coefficients
  1166        REAL(kind=dp), DIMENSION(3)                        :: ra
  1167        REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
  1168        TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
  1169        TYPE(cdft_control_type), POINTER                   :: cdft_control
  1170        TYPE(cell_type), POINTER                           :: cell

===== qs_cdft_methods.F:1230-1270 =====
  1230           CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
  1231           CALL rs_grid_zero(rs_rho_constr)
  1232           ! Compute Gaussian density over single atoms (rs_single) when atomic charges are requested
  1233           IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
  1234              ALLOCATE (rs_single(cdft_control%natoms))
  1235              DO i = 1, cdft_control%natoms
  1236                 NULLIFY (rs_single(i)%rs_grid)
  1237                 CALL rs_grid_create(rs_single(i)%rs_grid, auxbas_rs_desc)
  1238                 CALL rs_grid_zero(rs_single(i)%rs_grid)
  1239                 compute_charge(cdft_control%atoms(i)) = .TRUE.
  1240              END DO
  1241           END IF
  1242  
  1243           ! Collocate Gaussians
  1244           eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
  1245           ALLOCATE (pab(1, 1))
  1246           nthread = 1
  1247           ithread = 0
  1248  
  1249           DO ikind = 1, SIZE(atomic_kind_set)
  1250              numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
  1251              IF (numexp <= 0) CYCLE
  1252              CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
  1253              ALLOCATE (cores(natom))
  1254  
  1255              DO iex = 1, numexp
  1256                 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
  1257                 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
  1258                 npme = 0
  1259                 cores = 0
  1260                 DO iatom = 1, natom
  1261                    atom_a = atom_list(iatom)
  1262                    ra(:) = pbc(particle_set(atom_a)%r, cell)
  1263                    IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
  1264                       ! replicated realspace grid, split the atoms up between procs
  1265                       IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
  1266                          npme = npme + 1
  1267                          cores(npme) = iatom
  1268                       ENDIF
  1269                    ELSE
  1270                       npme = npme + 1

===== qs_cdft_methods.F:1320-1365 =====
  1320              DEALLOCATE (cores)
  1321           END DO
  1322           DEALLOCATE (pab)
  1323  
  1324           ! Transfer rs_rho_all to the correct grid and save it
  1325           CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
  1326           IF (igroup == 1) THEN
  1327              IF (ASSOCIATED(fnorm)) THEN
  1328                 CALL pw_pool_give_back_pw(auxbas_pw_pool, fnorm%pw)
  1329              END IF
  1330              ALLOCATE (fnorm)
  1331              CALL pw_pool_create_pw(auxbas_pw_pool, fnorm%pw, use_data=REALDATA3D, &
  1332                                     in_space=REALSPACE)
  1333              CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
  1334              CALL rs_pw_transfer(rs_rho_all, fnorm%pw, rs2pw)
  1335              CALL rs_grid_release(rs_rho_all)
  1336           END IF
  1337           ! Compute CDFT weight function
  1338           CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, &
  1339                                  in_space=REALSPACE)
  1340           CALL rs_pw_transfer(rs_rho_constr, tmp%pw, rs2pw)
  1341           CALL rs_grid_release(rs_rho_constr)
  1342           CALL hfun_scale(cdft_control%group(igroup)%weight%pw%cr3d, tmp%pw%cr3d, &
  1343                           fnorm%pw%cr3d, divide=.TRUE.)
  1344           CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw)
  1345           ! Compute atomic weight functions if charges are needed
  1346           IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
  1347              CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, &
  1348                                     in_space=REALSPACE)
  1349              DO i = 1, cdft_control%natoms
  1350                 CALL rs_pw_transfer(rs_single(i)%rs_grid, tmp%pw, rs2pw)
  1351                 CALL rs_grid_release(rs_single(i)%rs_grid)
  1352                 CALL hfun_scale(cdft_control%charge(i)%pw%cr3d, tmp%pw%cr3d, &
  1353                                 fnorm%pw%cr3d, divide=.TRUE.)
  1354              END DO
  1355              CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw)
  1356              DEALLOCATE (rs_single)
  1357              DEALLOCATE (compute_charge)
  1358           END IF
  1359        END DO
  1360  
  1361        DEALLOCATE (is_constraint)
  1362        DEALLOCATE (coefficients)
  1363        CALL timestop(handle)
  1364  
  1365     END SUBROUTINE hirshfeld_constraint_low

===== qs_ks_apply_restraints.F:62-128 =====
    62     SUBROUTINE qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
    63        TYPE(qs_environment_type), POINTER                 :: qs_env
    64        TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    65        LOGICAL, INTENT(in)                                :: calculate_forces
    66        TYPE(cdft_control_type), POINTER                   :: cdft_control
    67  
    68        INTEGER                                            :: iatom, igroup, natom
    69        LOGICAL                                            :: do_kpoints
    70        REAL(KIND=dp)                                      :: inv_vol
    71        TYPE(dft_control_type), POINTER                    :: dft_control
    72  
    73        NULLIFY (dft_control)
    74        CALL get_qs_env(qs_env, dft_control=dft_control)
    75        IF (dft_control%qs_control%cdft) THEN
    76           cdft_control => dft_control%qs_control%cdft_control
    77           ! Test no k-points
    78           CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
    79           IF (do_kpoints) CPABORT("CDFT constraints with k-points not supported.")
    80  
    81           SELECT CASE (cdft_control%type)
    82           CASE (outer_scf_becke_constraint, outer_scf_hirshfeld_constraint)
    83              IF (cdft_control%need_pot) THEN
    84                 ! First SCF iteraration => allocate storage
    85                 DO igroup = 1, SIZE(cdft_control%group)
    86                    CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%weight%pw, &
    87                                           use_data=REALDATA3D, in_space=REALSPACE)
    88                    ! Sanity check
    89                    IF (cdft_control%group(igroup)%constraint_type /= cdft_charge_constraint &
    90                        .AND. dft_control%nspins == 1) &
    91                       CALL cp_abort(__LOCATION__, &
    92                                     "Spin constraints require a spin polarized calculation.")
    93                 END DO
    94                 IF (cdft_control%atomic_charges) THEN
    95                    IF (.NOT. ASSOCIATED(cdft_control%charge)) &
    96                       ALLOCATE (cdft_control%charge(cdft_control%natoms))
    97                    DO iatom = 1, cdft_control%natoms
    98                       CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%charge(iatom)%pw, &
    99                                              use_data=REALDATA3D, in_space=REALSPACE)
   100                    END DO
   101                 END IF
   102                 ! Another sanity check
   103                 CALL get_qs_env(qs_env, natom=natom)
   104                 IF (natom < cdft_control%natoms) &
   105                    CALL cp_abort(__LOCATION__, &
   106                                  "The number of constraint atoms exceeds the total number of atoms.")
   107              ELSE
   108                 DO igroup = 1, SIZE(cdft_control%group)
   109                    inv_vol = 1.0_dp/cdft_control%group(igroup)%weight%pw%pw_grid%dvol
   110                    CALL pw_scale(cdft_control%group(igroup)%weight%pw, inv_vol)
   111                 END DO
   112              END IF
   113              ! Build/Integrate CDFT constraints with selected population analysis method
   114              IF (cdft_control%type == outer_scf_becke_constraint) THEN
   115                 CALL becke_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
   116              ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
   117                 CALL hirshfeld_constraint(qs_env, calc_pot=cdft_control%need_pot, calculate_forces=calculate_forces)
   118              END IF
   119              DO igroup = 1, SIZE(cdft_control%group)
   120                 CALL pw_scale(cdft_control%group(igroup)%weight%pw, cdft_control%group(igroup)%weight%pw%pw_grid%dvol)
   121              END DO
   122              IF (cdft_control%need_pot) cdft_control%need_pot = .FALSE.
   123           CASE DEFAULT
   124              CPABORT("Unknown constraint type.")
   125           END SELECT
   126        END IF
   127  
   128     END SUBROUTINE qs_ks_cdft_constraint

===== qs_scf.F:992-1055 =====
   992        cdft_outer_loop: DO
   993           ! Change outer_scf settings to OT settings
   994           CALL outer_loop_switch(scf_env, scf_control, cdft_control, cdft2ot)
   995           ! Solve electronic structure with fixed value of constraint
   996           CALL scf_env_do_scf(scf_env=scf_env, scf_control=scf_control, qs_env=qs_env, &
   997                               converged=converged, should_stop=should_stop)
   998           ! Decide whether to reuse the preconditioner on the next iteration
   999           IF (cdft_control%reuse_precond) THEN
  1000              ! For convergence in exactly one step, the preconditioner is always reused (assuming max_reuse > 0)
  1001              ! usually this means that the electronic structure has already converged to the correct state
  1002              ! but the constraint optimizer keeps jumping over the optimal solution
  1003              IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count == 1 &
  1004                  .AND. cdft_control%total_steps /= 1) &
  1005                 cdft_control%nreused = cdft_control%nreused - 1
  1006              ! SCF converged in less than precond_freq steps
  1007              IF (scf_env%outer_scf%iter_count == 1 .AND. scf_env%iter_count .LE. cdft_control%precond_freq .AND. &
  1008                  cdft_control%total_steps /= 1 .AND. cdft_control%nreused .LT. cdft_control%max_reuse) THEN
  1009                 reuse_precond = .TRUE.
  1010                 cdft_control%nreused = cdft_control%nreused + 1
  1011              ELSE
  1012                 reuse_precond = .FALSE.
  1013                 cdft_control%nreused = 0
  1014              END IF
  1015           END IF
  1016           ! Update history purging counters
  1017           IF (first_iteration .AND. cdft_control%purge_history) THEN
  1018              cdft_control%istep = cdft_control%istep + 1
  1019              IF (scf_env%outer_scf%iter_count .GT. 1) THEN
  1020                 cdft_control%nbad_conv = cdft_control%nbad_conv + 1
  1021                 IF (cdft_control%nbad_conv .GE. cdft_control%purge_freq .AND. &
  1022                     cdft_control%istep .GE. cdft_control%purge_offset) THEN
  1023                    cdft_control%nbad_conv = 0
  1024                    cdft_control%istep = 0
  1025                    cdft_control%should_purge = .TRUE.
  1026                 END IF
  1027              END IF
  1028           END IF
  1029           first_iteration = .FALSE.
  1030           ! Change outer_scf settings to CDFT settings
  1031           CALL outer_loop_switch(scf_env, scf_control, cdft_control, ot2cdft)
  1032           CALL qs_scf_check_outer_exit(qs_env, scf_env, scf_control, should_stop, &
  1033                                        cdft_loop_converged, exit_cdft_loop)
  1034           CALL qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
  1035                                 energy, cdft_control%total_steps, &
  1036                                 should_stop, cdft_loop_converged, cdft_loop=.TRUE.)
  1037           IF (exit_cdft_loop) EXIT cdft_outer_loop
  1038           ! Check if the inverse Jacobian needs to be calculated
  1039           CALL qs_calculate_inverse_jacobian(qs_env)
  1040           ! Check if a line search should be performed to find an optimal step size for the optimizer
  1041           CALL qs_cdft_line_search(qs_env)
  1042           ! Optimize constraint
  1043           CALL outer_loop_optimize(scf_env, scf_control)
  1044           CALL outer_loop_update_qs_env(qs_env, scf_env)
  1045           CALL qs_ks_did_change(ks_env, potential_changed=.TRUE.)
  1046        END DO cdft_outer_loop
  1047  
  1048        cdft_control%ienergy = cdft_control%ienergy + 1
  1049  
  1050        ! Store needed arrays for ET coupling calculation
  1051        IF (cdft_control%do_et) THEN
  1052           CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, mos=mos)
  1053           nvar = SIZE(cdft_control%target)
  1054           ! Matrix representation of weight function
  1055           ALLOCATE (cdft_control%wmat(nvar))

===== qs_ks_methods.F:470-490 =====
   470              energy%hartree = energy%hartree + energy%surf_dipole
   471           END IF
   472        END IF
   473  
   474        ! SIC
   475        CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, &
   476                               just_energy, calculate_forces, auxbas_pw_pool)
   477  
   478        IF (gapw) THEN
   479           CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
   480           CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
   481        END IF
   482  
   483        ! Check if CDFT constraint is needed
   484        CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
   485  
   486        ! Adds the External Potential if requested
   487        IF (dft_control%apply_external_potential) THEN
   488           ! Compute the energy due to the external potential
   489           ee_ener = 0.0_dp
   490           DO ispin = 1, nspins

===== qs_ks_utils.F:1321-1345 =====
  1321              ! CDFT constraint contribution
  1322              IF (dft_control%qs_control%cdft) THEN
  1323                 DO igroup = 1, SIZE(cdft_control%group)
  1324                    SELECT CASE (cdft_control%group(igroup)%constraint_type)
  1325                    CASE (cdft_charge_constraint)
  1326                       sign = 1.0_dp
  1327                    CASE (cdft_magnetization_constraint)
  1328                       IF (ispin == 1) THEN
  1329                          sign = 1.0_dp
  1330                       ELSE
  1331                          sign = -1.0_dp
  1332                       END IF
  1333                    CASE (cdft_alpha_constraint)
  1334                       sign = 1.0_dp
  1335                       IF (ispin == 2) CYCLE
  1336                    CASE (cdft_beta_constraint)
  1337                       sign = 1.0_dp
  1338                       IF (ispin == 1) CYCLE
  1339                    CASE DEFAULT
  1340                       CPABORT("Unknown constraint type.")
  1341                    END SELECT
  1342                    v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d &
  1343                                                  +sign*cdft_control%group(igroup)%weight%pw%cr3d* &
  1344                                                  cdft_control%strength(igroup)
  1345                 END DO

===== qs_outer_scf.F:102-175 =====
   102     SUBROUTINE outer_loop_gradient(qs_env, scf_env)
   103        TYPE(qs_environment_type), POINTER                 :: qs_env
   104        TYPE(qs_scf_env_type), POINTER                     :: scf_env
   105  
   106        CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_gradient'
   107  
   108        INTEGER                                            :: handle, ihistory, ivar, n
   109        LOGICAL                                            :: is_constraint
   110        TYPE(cdft_control_type), POINTER                   :: cdft_control
   111        TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
   112        TYPE(dft_control_type), POINTER                    :: dft_control
   113        TYPE(qs_energy_type), POINTER                      :: energy
   114        TYPE(s2_restraint_type), POINTER                   :: s2_restraint_control
   115        TYPE(scf_control_type), POINTER                    :: scf_control
   116  
   117        CALL timeset(routineN, handle)
   118  
   119        CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, &
   120                        dft_control=dft_control, energy=energy)
   121        CPASSERT(scf_control%outer_scf%have_scf)
   122  
   123        ihistory = scf_env%outer_scf%iter_count
   124        CPASSERT(ihistory <= SIZE(scf_env%outer_scf%energy, 1))
   125  
   126        scf_env%outer_scf%energy(ihistory) = energy%total
   127  
   128        SELECT CASE (scf_control%outer_scf%type)
   129        CASE (outer_scf_none)
   130           ! just pass the inner loop scf criterion to the outer loop one
   131           scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta
   132           scf_env%outer_scf%gradient(1, ihistory) = scf_env%iter_delta
   133        CASE (outer_scf_ddapc_constraint)
   134           CPASSERT(dft_control%qs_control%ddapc_restraint)
   135           DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
   136              NULLIFY (ddapc_restraint_control)
   137              ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control
   138              is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
   139              IF (is_constraint) EXIT
   140           END DO
   141           CPASSERT(is_constraint)
   142  
   143           scf_env%outer_scf%variables(:, ihistory) = ddapc_restraint_control%strength
   144           scf_env%outer_scf%gradient(:, ihistory) = ddapc_restraint_control%ddapc_order_p - &
   145                                                     ddapc_restraint_control%target
   146        CASE (outer_scf_s2_constraint)
   147           CPASSERT(dft_control%qs_control%s2_restraint)
   148           s2_restraint_control => dft_control%qs_control%s2_restraint_control
   149           is_constraint = (s2_restraint_control%functional_form == do_s2_constraint)
   150           CPASSERT(is_constraint)
   151  
   152           scf_env%outer_scf%variables(:, ihistory) = s2_restraint_control%strength
   153           scf_env%outer_scf%gradient(:, ihistory) = s2_restraint_control%s2_order_p - &
   154                                                     s2_restraint_control%target
   155        CASE (outer_scf_cdft_constraint)
   156           CPASSERT(dft_control%qs_control%cdft)
   157           cdft_control => dft_control%qs_control%cdft_control
   158           DO ivar = 1, SIZE(scf_env%outer_scf%gradient, 1)
   159              scf_env%outer_scf%variables(ivar, ihistory) = cdft_control%strength(ivar)
   160              scf_env%outer_scf%gradient(ivar, ihistory) = cdft_control%value(ivar) - &
   161                                                           cdft_control%target(ivar)
   162           END DO
   163        CASE (outer_scf_basis_center_opt)
   164           CALL qs_basis_center_gradient(qs_env)
   165           scf_env%outer_scf%gradient(:, ihistory) = return_basis_center_gradient_norm(qs_env)
   166  
   167        CASE DEFAULT
   168           CPABORT("")
   169  
   170        END SELECT
   171  
   172        CALL timestop(handle)
   173  
   174     END SUBROUTINE outer_loop_gradient
   175  

===== qs_outer_scf.F:386-430 =====
   386     SUBROUTINE outer_loop_update_qs_env(qs_env, scf_env)
   387        TYPE(qs_environment_type), POINTER                 :: qs_env
   388        TYPE(qs_scf_env_type), POINTER                     :: scf_env
   389  
   390        CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_update_qs_env'
   391  
   392        INTEGER                                            :: handle, ihistory, n
   393        LOGICAL                                            :: is_constraint
   394        TYPE(cdft_control_type), POINTER                   :: cdft_control
   395        TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
   396        TYPE(dft_control_type), POINTER                    :: dft_control
   397        TYPE(s2_restraint_type), POINTER                   :: s2_restraint_control
   398        TYPE(scf_control_type), POINTER                    :: scf_control
   399  
   400        CALL timeset(routineN, handle)
   401  
   402        CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
   403        ihistory = scf_env%outer_scf%iter_count
   404  
   405        SELECT CASE (scf_control%outer_scf%type)
   406        CASE (outer_scf_none)
   407           ! do nothing
   408        CASE (outer_scf_ddapc_constraint)
   409           DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
   410              NULLIFY (ddapc_restraint_control)
   411              ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control
   412              is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint)
   413              IF (is_constraint) EXIT
   414           END DO
   415           ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
   416        CASE (outer_scf_s2_constraint)
   417           s2_restraint_control => dft_control%qs_control%s2_restraint_control
   418           s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1)
   419        CASE (outer_scf_cdft_constraint)
   420           cdft_control => dft_control%qs_control%cdft_control
   421           cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1)
   422        CASE (outer_scf_basis_center_opt)
   423           CALL qs_update_basis_center_pos(qs_env)
   424        CASE DEFAULT
   425           CPABORT("")
   426        END SELECT
   427  
   428        CALL timestop(handle)
   429  
   430     END SUBROUTINE outer_loop_update_qs_env

===== hirshfeld_methods.F:141-224 =====
   141     SUBROUTINE create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, radius, radii_list)
   142        TYPE(hirshfeld_type), POINTER                      :: hirshfeld_env
   143        TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
   144        TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
   145        REAL(KIND=dp), OPTIONAL                            :: radius
   146        REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: radii_list
   147  
   148        INTEGER, PARAMETER                                 :: ngto = 8
   149  
   150        CHARACTER(len=2)                                   :: esym
   151        INTEGER                                            :: ikind, nkind
   152        LOGICAL                                            :: found
   153        REAL(KIND=dp)                                      :: al, rco, zeff
   154        REAL(KIND=dp), DIMENSION(ngto, 2)                  :: ppdens
   155        TYPE(atomic_kind_type), POINTER                    :: atomic_kind
   156        TYPE(qs_kind_type), POINTER                        :: qs_kind
   157  
   158        CPASSERT(ASSOCIATED(hirshfeld_env))
   159  
   160        nkind = SIZE(qs_kind_set)
   161        ALLOCATE (hirshfeld_env%kind_shape_fn(nkind))
   162  
   163        SELECT CASE (hirshfeld_env%shape_function_type)
   164        CASE (shape_function_gaussian)
   165           DO ikind = 1, nkind
   166              hirshfeld_env%kind_shape_fn(ikind)%numexp = 1
   167              ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(1))
   168              ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(1))
   169              CALL get_qs_kind(qs_kind_set(ikind), element_symbol=esym)
   170              rco = 2.0_dp
   171              SELECT CASE (hirshfeld_env%radius_type)
   172              CASE (radius_default)
   173                 CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found)
   174                 rco = MAX(rco, 1.0_dp)
   175              CASE (radius_user)
   176                 CPASSERT(PRESENT(radii_list))
   177                 CPASSERT(ASSOCIATED(radii_list))
   178                 CPASSERT(SIZE(radii_list) == nkind)
   179                 ! Note we assume that radii_list is correctly ordered
   180                 rco = radii_list(ikind)
   181              CASE (radius_vdw)
   182                 CALL get_ptable_info(symbol=esym, vdw_radius=rco, found=found)
   183                 IF (.NOT. found) THEN
   184                    rco = MAX(rco, 1.0_dp)
   185                 ELSE
   186                    IF (hirshfeld_env%use_bohr) &
   187                       rco = cp_unit_to_cp2k(rco, "angstrom")
   188                 END IF
   189              CASE (radius_covalent)
   190                 CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found)
   191                 IF (.NOT. found) THEN
   192                    rco = MAX(rco, 1.0_dp)
   193                 ELSE
   194                    IF (hirshfeld_env%use_bohr) &
   195                       rco = cp_unit_to_cp2k(rco, "angstrom")
   196                 END IF
   197              CASE (radius_single)
   198                 CPASSERT(PRESENT(radius))
   199                 rco = radius
   200              END SELECT
   201              al = 0.5_dp/rco**2
   202              hirshfeld_env%kind_shape_fn(ikind)%zet(1) = al
   203              hirshfeld_env%kind_shape_fn(ikind)%coef(1) = (al/pi)**1.5_dp
   204           END DO
   205        CASE (shape_function_density)
   206           ! calculate atomic density
   207           DO ikind = 1, nkind
   208              atomic_kind => atomic_kind_set(ikind)
   209              qs_kind => qs_kind_set(ikind)
   210              CALL calculate_atomic_density(ppdens(:, :), atomic_kind, qs_kind, ngto, &
   211                                            confine=.FALSE.)
   212              hirshfeld_env%kind_shape_fn(ikind)%numexp = ngto
   213              ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(ngto))
   214              ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(ngto))
   215              hirshfeld_env%kind_shape_fn(ikind)%zet(:) = ppdens(:, 1)
   216              CALL get_qs_kind(qs_kind, zeff=zeff)
   217              hirshfeld_env%kind_shape_fn(ikind)%coef(:) = ppdens(:, 2)/zeff
   218           END DO
   219  
   220        CASE DEFAULT
   221           CPABORT("Unknown shape function")
   222        END SELECT
   223  
   224     END SUBROUTINE create_shape_function
