comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Bump fortitude to version 0.9.0 and consider some of its new rules (#5597)

Browse source
This commit is contained in:
Frederick Stein 2026-07-17 08:19:44 +02:00 • committed by GitHub
parent 39fe4c47c6
commit 3caaa28e3b
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
622 changed files with 7366 additions and 7363 deletions
Download patch file Download diff file Expand all files Collapse all files

15
src/accint_weights_forces.F
View file

@ -412,7 +412,7 @@ CONTAINS
NULLIFY (rho_r, rho_g, tau_r, tau_g)
IF (rho_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, rho_g_base, rho_r, rho_g)
ELSEIF (ASSOCIATED(rho_r_base)) THEN
ELSE IF (ASSOCIATED(rho_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_r_base, rho_r, rho_g)
ELSE
CPABORT("Fine Grid in xc_density requires rho_r or rho_g")
@ -420,7 +420,7 @@ CONTAINS
IF (rho_tau_valid) THEN
IF (rho_tau_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, tau_g_base, tau_r, tau_g)
ELSEIF (ASSOCIATED(tau_r_base)) THEN
ELSE IF (ASSOCIATED(tau_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_r_base, tau_r, tau_g)
ELSE
CPABORT("Fine Grid in xc_density requires tau_r or tau_g")
@ -482,7 +482,7 @@ CONTAINS
NULLIFY (rho1_r, rho1_g, tau1_r, tau1_g)
IF (rho1_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, rho1_g_base, rho1_r, rho1_g)
ELSEIF (ASSOCIATED(rho1_r_base)) THEN
ELSE IF (ASSOCIATED(rho1_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho1_r_base, rho1_r, rho1_g)
ELSE
CPABORT("Fine Grid in xc_density requires rho1_r or rho1_g")
@ -490,7 +490,7 @@ CONTAINS
IF (rho1_tau_valid) THEN
IF (rho1_tau_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, tau1_g_base, tau1_r, tau1_g)
ELSEIF (ASSOCIATED(tau1_r_base)) THEN
ELSE IF (ASSOCIATED(tau1_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau1_r_base, tau1_r, tau1_g)
ELSE
CPABORT("Fine Grid in xc_density requires tau1_r or tau1_g")
@ -518,7 +518,7 @@ CONTAINS
NULLIFY (rho1_r, rho1_g, tau1_r, tau1_g)
IF (rho1_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, rho1_g_base, rho1_r, rho1_g)
ELSEIF (ASSOCIATED(rho1_r_base)) THEN
ELSE IF (ASSOCIATED(rho1_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho1_r_base, rho1_r, rho1_g)
ELSE
CPABORT("Fine Grid in xc_density requires rho1_r or rho1_g")
@ -526,7 +526,7 @@ CONTAINS
IF (rho1_tau_valid) THEN
IF (rho1_tau_g_valid) THEN
CALL create_density_on_pool(xc_pw_pool, tau1_g_base, tau1_r, tau1_g)
ELSEIF (ASSOCIATED(tau1_r_base)) THEN
ELSE IF (ASSOCIATED(tau1_r_base)) THEN
CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau1_r_base, tau1_r, tau1_g)
ELSE
CPABORT("Fine Grid in xc_density requires tau1_r or tau1_g")
@ -574,8 +574,9 @@ CONTAINS
DEALLOCATE (vxc_rho)
END IF
IF (ASSOCIATED(vxc_tau)) THEN
IF (.NOT. ASSOCIATED(tau1_r)) &
IF (.NOT. ASSOCIATED(tau1_r)) THEN
CPABORT("Tau response density required for mGGA xc_density")
END IF
DO ispin = 1, nspins
CALL pw_multiply_with(vxc_tau(ispin), tau1_r(ispin))
CALL pw_axpy(vxc_tau(ispin), exc, 1.0_dp)

12
src/admm_dm_methods.F
View file

@ -76,8 +76,9 @@ CONTAINS
CPABORT("admm_dm_calc_rho_aux: unknown method")
END SELECT
IF (admm_dm%purify) &
IF (admm_dm%purify) THEN
CALL purify_mcweeny(qs_env)
END IF
CALL update_rho_aux(qs_env)
@ -120,8 +121,9 @@ CONTAINS
CPABORT("admm_dm_merge_ks_matrix: unknown method")
END SELECT
IF (admm_dm%purify) &
IF (admm_dm%purify) THEN
CALL dbcsr_deallocate_matrix_set(matrix_ks_merge)
END IF
CALL timestop(handle)
@ -217,8 +219,9 @@ CONTAINS
IF (admm_dm%block_map(iatom, jatom) == 1) THEN
CALL dbcsr_get_block_p(rho_ao_aux(ispin)%matrix, &
row=iatom, col=jatom, BLOCK=sparse_block_aux, found=found)
IF (found) &
IF (found) THEN
sparse_block_aux = sparse_block
END IF
END IF
END DO
CALL dbcsr_iterator_stop(iter)
@ -333,8 +336,9 @@ CONTAINS
CALL dbcsr_iterator_start(iter, matrix_ks_merge(ispin)%matrix)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, iatom, jatom, sparse_block)
IF (admm_dm%block_map(iatom, jatom) == 0) &
IF (admm_dm%block_map(iatom, jatom) == 0) THEN
sparse_block = 0.0_dp
END IF
END DO
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_ks_merge(ispin)%matrix, 1.0_dp, 1.0_dp)

3
src/admm_dm_types.F
View file

@ -105,8 +105,9 @@ CONTAINS
DEALLOCATE (admm_dm%matrix_a)
END IF
IF (ASSOCIATED(admm_dm%block_map)) &
IF (ASSOCIATED(admm_dm%block_map)) THEN
DEALLOCATE (admm_dm%block_map)
END IF
DEALLOCATE (admm_dm%mcweeny_history)
DEALLOCATE (admm_dm)

15
src/admm_methods.F
View file

@ -226,12 +226,13 @@ CONTAINS
END IF
IF (admm_env%purification_method == do_admm_purify_cauchy) &
IF (admm_env%purification_method == do_admm_purify_cauchy) THEN
CALL purify_dm_cauchy(admm_env, &
mo_set=mos_aux_fit(ispin), &
density_matrix=rho_ao_aux(ispin)%matrix, &
ispin=ispin, &
blocked=admm_env%block_dm)
END IF
!GPW is the default, PW density is computed using the AUX_FIT basis and task_list
!If GAPW, the we use the AUX_FIT_SOFT basis and task list
@ -2180,13 +2181,15 @@ CONTAINS
IF (my_kpgrp) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
IF (.NOT. use_real_wfn) &
IF (.NOT. use_real_wfn) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, &
para_env, info(indx, 2))
END IF
ELSE
CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
IF (.NOT. use_real_wfn) &
IF (.NOT. use_real_wfn) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, fmdummy, para_env, info(indx, 2))
END IF
END IF
END DO
END DO
@ -2984,12 +2987,14 @@ CONTAINS
IF (my_kpgrp) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux, work_aux_aux, para_env, info(indx, 1))
IF (.NOT. use_real_wfn) &
IF (.NOT. use_real_wfn) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, work_aux_aux2, para_env, info(indx, 2))
END IF
ELSE
CALL cp_fm_start_copy_general(admm_env%work_aux_aux, fmdummy, para_env, info(indx, 1))
IF (.NOT. use_real_wfn) &
IF (.NOT. use_real_wfn) THEN
CALL cp_fm_start_copy_general(admm_env%work_aux_aux2, fmdummy, para_env, info(indx, 2))
END IF
END IF
END DO
END DO

18
src/admm_types.F
View file

@ -382,12 +382,15 @@ CONTAINS
admm_env%aux_x_param(:) = admm_control%aux_x_param(:)
!ADMMP, ADMMQ, ADMMS
IF ((.NOT. admm_env%charge_constrain) .AND. (admm_env%scaling_model == do_admm_exch_scaling_merlot)) &
IF ((.NOT. admm_env%charge_constrain) .AND. (admm_env%scaling_model == do_admm_exch_scaling_merlot)) THEN
admm_env%do_admmp = .TRUE.
IF (admm_env%charge_constrain .AND. (admm_env%scaling_model == do_admm_exch_scaling_none)) &
END IF
IF (admm_env%charge_constrain .AND. (admm_env%scaling_model == do_admm_exch_scaling_none)) THEN
admm_env%do_admmq = .TRUE.
IF (admm_env%charge_constrain .AND. (admm_env%scaling_model == do_admm_exch_scaling_merlot)) &
END IF
IF (admm_env%charge_constrain .AND. (admm_env%scaling_model == do_admm_exch_scaling_merlot)) THEN
admm_env%do_admms = .TRUE.
END IF
IF ((admm_control%method == do_admm_blocking_purify_full) .OR. &
(admm_control%method == do_admm_blocked_projection)) THEN
@ -483,13 +486,16 @@ CONTAINS
DEALLOCATE (admm_env%eigvals_lambda)
DEALLOCATE (admm_env%eigvals_P_to_be_purified)
IF (ASSOCIATED(admm_env%block_map)) &
IF (ASSOCIATED(admm_env%block_map)) THEN
DEALLOCATE (admm_env%block_map)
END IF
IF (ASSOCIATED(admm_env%xc_section_primary)) &
IF (ASSOCIATED(admm_env%xc_section_primary)) THEN
CALL section_vals_release(admm_env%xc_section_primary)
IF (ASSOCIATED(admm_env%xc_section_aux)) &
END IF
IF (ASSOCIATED(admm_env%xc_section_aux)) THEN
CALL section_vals_release(admm_env%xc_section_aux)
END IF
IF (ASSOCIATED(admm_env%admm_gapw_env)) CALL admm_gapw_env_release(admm_env%admm_gapw_env)
IF (ASSOCIATED(admm_env%admm_dm)) CALL admm_dm_release(admm_env%admm_dm)

35
src/almo_scf.F
View file

@ -1400,7 +1400,6 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_delocalization'
INTEGER :: handle, ispin, unit_nr
LOGICAL :: almo_experimental
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: no_quench
TYPE(optimizer_options_type) :: arbitrary_optimizer
@ -1462,25 +1461,6 @@ CONTAINS
!!!! are commented out at the moment because some of their
!!!! routines have not been thoroughly tested.
!!!! if we have virtuals pre-optimize and truncate them
!!!IF (almo_scf_env%need_virtuals) THEN
!!! SELECT CASE (almo_scf_env%deloc_truncate_virt)
!!! CASE (virt_full)
!!! ! simply copy virtual orbitals from matrix_v_full_blk to matrix_v_blk
!!! DO ispin=1,almo_scf_env%nspins
!!! CALL dbcsr_copy(almo_scf_env%matrix_v_blk(ispin),&
!!! almo_scf_env%matrix_v_full_blk(ispin))
!!! ENDDO
!!! CASE (virt_number,virt_occ_size)
!!! CALL split_v_blk(almo_scf_env)
!!! !CALL truncate_subspace_v_blk(qs_env,almo_scf_env)
!!! CASE DEFAULT
!!! CPErrorMessage(cp_failure_level,routineP,"illegal method for virtual space truncation")
!!! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
!!! END SELECT
!!!ENDIF
!!!CALL harris_foulkes_correction(qs_env,almo_scf_env)
IF (almo_scf_env%xalmo_update_algorithm == almo_scf_pcg) THEN
CALL almo_scf_xalmo_pcg(qs_env=qs_env, &
@ -1585,18 +1565,6 @@ CONTAINS
perturbation_only=.FALSE., &
special_case=xalmo_case_normal)
! RZK-warning THIS IS A HACK TO GET ORBITAL ENERGIES
almo_experimental = .FALSE.
IF (almo_experimental) THEN
almo_scf_env%perturbative_delocalization = .TRUE.
!DO ispin=1,almo_scf_env%nspins
! CALL dbcsr_copy(almo_scf_env%matrix_t(ispin),&
! almo_scf_env%matrix_t_blk(ispin))
!ENDDO
CALL almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, &
arbitrary_optimizer)
END IF ! experimental
ELSE IF (almo_scf_env%xalmo_update_algorithm == almo_scf_trustr) THEN
CALL almo_scf_xalmo_trustr(qs_env=qs_env, &
@ -2622,8 +2590,9 @@ CONTAINS
ALLOCATE (almo_scf_env%matrix_p_blk(nspins))
ALLOCATE (almo_scf_env%matrix_ks(nspins))
ALLOCATE (almo_scf_env%matrix_ks_blk(nspins))
IF (almo_scf_env%need_previous_ks) &
IF (almo_scf_env%need_previous_ks) THEN
ALLOCATE (almo_scf_env%matrix_ks_0deloc(nspins))
END IF
DO ispin = 1, nspins
! RZK-warning copy with symmery but remember that this might cause problems
CALL dbcsr_create(almo_scf_env%matrix_p(ispin), &

16
src/almo_scf_diis_types.F
View file

@ -258,8 +258,9 @@ CONTAINS
! update the buffer length
old_buffer_length = diis_env%buffer_length
diis_env%buffer_length = diis_env%buffer_length + 1
IF (diis_env%buffer_length > diis_env%max_buffer_length) &
IF (diis_env%buffer_length > diis_env%max_buffer_length) THEN
diis_env%buffer_length = diis_env%max_buffer_length
END IF
!!!! resize B matrix
!!!IF (old_buffer_length.lt.diis_env%buffer_length) THEN
@ -402,24 +403,12 @@ CONTAINS
! use the eigensystem to invert (implicitly) B matrix
! and compute the extrapolation coefficients
!! ALLOCATE(tmp1(diis_env%buffer_length+1,1))
!! ALLOCATE(coeff(diis_env%buffer_length+1,1))
!! tmp1(:,1)=-1.0_dp*m_b_copy(1,:)/eigenvalues(:)
!! coeff=MATMUL(m_b_copy,tmp1)
!! DEALLOCATE(tmp1)
ALLOCATE (tmp1(diis_env%buffer_length + 1))
ALLOCATE (coeff(diis_env%buffer_length + 1))
tmp1(:) = -1.0_dp*m_b_copy(1, :)/eigenvalues(:)
coeff(:) = MATMUL(m_b_copy, tmp1)
DEALLOCATE (tmp1)
!IF (unit_nr.gt.0) THEN
! DO im=1,diis_env%buffer_length+1
! WRITE(unit_nr,*) diis_env%m_b(idomain)%mdata(im,:)
! ENDDO
! WRITE (unit_nr,*) coeff(:,1)
!ENDIF
! extrapolate the variable
checksum = 0.0_dp
IF (diis_env%diis_env_type == diis_env_dbcsr) THEN
@ -441,7 +430,6 @@ CONTAINS
checksum = checksum + coeff(im + 1)
END DO
END IF
!WRITE(*,*) checksum
DEALLOCATE (coeff)

102
src/almo_scf_env_methods.F
View file

@ -364,108 +364,6 @@ CONTAINS
almo_scf_env%activate = 0
END IF
!CALL section_vals_val_get(almo_scf_section,"DOMAIN_LAYOUT_AOS",&
! i_val=almo_scf_env%domain_layout_aos)
!CALL section_vals_val_get(almo_scf_section,"DOMAIN_LAYOUT_MOS",&
! i_val=almo_scf_env%domain_layout_mos)
!CALL section_vals_val_get(almo_scf_section,"MATRIX_CLUSTERING_AOS",&
! i_val=almo_scf_env%mat_distr_aos)
!CALL section_vals_val_get(almo_scf_section,"MATRIX_CLUSTERING_MOS",&
! i_val=almo_scf_env%mat_distr_mos)
!CALL section_vals_val_get(almo_scf_section,"CONSTRAINT_TYPE",&
! i_val=almo_scf_env%constraint_type)
!CALL section_vals_val_get(almo_scf_section,"MU",&
! r_val=almo_scf_env%mu)
!CALL section_vals_val_get(almo_scf_section,"FIXED_MU",&
! l_val=almo_scf_env%fixed_mu)
!CALL section_vals_val_get(almo_scf_section,"EPS_USE_PREV_AS_GUESS",&
! r_val=almo_scf_env%eps_prev_guess)
!CALL section_vals_val_get(almo_scf_section,"MIXING_FRACTION",&
! r_val=almo_scf_env%mixing_fraction)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_TENSOR_TYPE",&
! i_val=almo_scf_env%deloc_cayley_tensor_type)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_CONJUGATOR",&
! i_val=almo_scf_env%deloc_cayley_conjugator)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_MAX_ITER",&
! i_val=almo_scf_env%deloc_cayley_max_iter)
!CALL section_vals_val_get(almo_scf_section,"DELOC_USE_OCC_ORBS",&
! l_val=almo_scf_env%deloc_use_occ_orbs)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_USE_VIRT_ORBS",&
! l_val=almo_scf_env%deloc_cayley_use_virt_orbs)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_LINEAR",&
! l_val=almo_scf_env%deloc_cayley_linear)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_EPS_CONVERGENCE",&
! r_val=almo_scf_env%deloc_cayley_eps_convergence)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_OCC_PRECOND",&
! l_val=almo_scf_env%deloc_cayley_occ_precond)
!CALL section_vals_val_get(almo_scf_section,"DELOC_CAYLEY_VIR_PRECOND",&
! l_val=almo_scf_env%deloc_cayley_vir_precond)
!CALL section_vals_val_get(almo_scf_section,"ALMO_UPDATE_ALGORITHM_BD",&
! i_val=almo_scf_env%almo_update_algorithm)
!CALL section_vals_val_get(almo_scf_section,"DELOC_TRUNCATE_VIRTUALS",&
! i_val=almo_scf_env%deloc_truncate_virt)
!CALL section_vals_val_get(almo_scf_section,"DELOC_VIRT_PER_DOMAIN",&
! i_val=almo_scf_env%deloc_virt_per_domain)
!
!CALL section_vals_val_get(almo_scf_section,"OPT_K_EPS_CONVERGENCE",&
! r_val=almo_scf_env%opt_k_eps_convergence)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_MAX_ITER",&
! i_val=almo_scf_env%opt_k_max_iter)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_OUTER_MAX_ITER",&
! i_val=almo_scf_env%opt_k_outer_max_iter)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_TRIAL_STEP_SIZE",&
! r_val=almo_scf_env%opt_k_trial_step_size)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_CONJUGATOR",&
! i_val=almo_scf_env%opt_k_conjugator)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_TRIAL_STEP_SIZE_MULTIPLIER",&
! r_val=almo_scf_env%opt_k_trial_step_size_multiplier)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_CONJ_ITER_START",&
! i_val=almo_scf_env%opt_k_conj_iter_start)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_PREC_ITER_START",&
! i_val=almo_scf_env%opt_k_prec_iter_start)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_CONJ_ITER_FREQ_RESET",&
! i_val=almo_scf_env%opt_k_conj_iter_freq)
!CALL section_vals_val_get(almo_scf_section,"OPT_K_PREC_ITER_FREQ_UPDATE",&
! i_val=almo_scf_env%opt_k_prec_iter_freq)
!
!CALL section_vals_val_get(almo_scf_section,"QUENCHER_RADIUS_TYPE",&
! i_val=almo_scf_env%quencher_radius_type)
!CALL section_vals_val_get(almo_scf_section,"QUENCHER_R0_FACTOR",&
! r_val=almo_scf_env%quencher_r0_factor)
!CALL section_vals_val_get(almo_scf_section,"QUENCHER_R1_FACTOR",&
! r_val=almo_scf_env%quencher_r1_factor)
!!CALL section_vals_val_get(almo_scf_section,"QUENCHER_R0_SHIFT",&
!! r_val=almo_scf_env%quencher_r0_shift)
!!
!!CALL section_vals_val_get(almo_scf_section,"QUENCHER_R1_SHIFT",&
!! r_val=almo_scf_env%quencher_r1_shift)
!!
!!almo_scf_env%quencher_r0_shift = cp_unit_to_cp2k(&
!! almo_scf_env%quencher_r0_shift,"angstrom")
!!almo_scf_env%quencher_r1_shift = cp_unit_to_cp2k(&
!! almo_scf_env%quencher_r1_shift,"angstrom")
!
!CALL section_vals_val_get(almo_scf_section,"QUENCHER_AO_OVERLAP_0",&
! r_val=almo_scf_env%quencher_s0)
!CALL section_vals_val_get(almo_scf_section,"QUENCHER_AO_OVERLAP_1",&
! r_val=almo_scf_env%quencher_s1)
!CALL section_vals_val_get(almo_scf_section,"ENVELOPE_AMPLITUDE",&
! r_val=almo_scf_env%envelope_amplitude)
!! how to read lists
!CALL section_vals_val_get(almo_scf_section,"INT_LIST01", &
! n_rep_val=n_rep)
!counter_i = 0
!DO k = 1,n_rep
! CALL section_vals_val_get(almo_scf_section,"INT_LIST01",&
! i_rep_val=k,i_vals=tmplist)
! DO jj = 1,SIZE(tmplist)
! counter_i=counter_i+1
! almo_scf_env%charge_of_domain(counter_i)=tmplist(jj)
! ENDDO
!ENDDO
almo_scf_env%domain_layout_aos = almo_domain_layout_molecular
almo_scf_env%domain_layout_mos = almo_domain_layout_molecular
almo_scf_env%mat_distr_aos = almo_mat_distr_molecular

93
src/almo_scf_methods.F
View file

@ -1620,8 +1620,9 @@ CONTAINS
my_algorithm = 0
IF (PRESENT(algorithm)) my_algorithm = algorithm
IF (my_algorithm == 1 .AND. (.NOT. PRESENT(para_env) .OR. .NOT. PRESENT(blacs_env))) &
IF (my_algorithm == 1 .AND. (.NOT. PRESENT(para_env) .OR. .NOT. PRESENT(blacs_env))) THEN
CPABORT("PARA and BLACS env are necessary for cholesky algorithm")
END IF
use_sigma_inv_guess = .FALSE.
IF (PRESENT(use_guess)) THEN
@ -2092,7 +2093,6 @@ CONTAINS
ALLOCATE (subm_in(ndomains))
ALLOCATE (subm_temp(ndomains))
ALLOCATE (subm_out(ndomains))
!!!TRIM ALLOCATE(subm_trimmer(ndomains))
CALL init_submatrices(subm_in)
CALL init_submatrices(subm_temp)
CALL init_submatrices(subm_out)
@ -2100,11 +2100,6 @@ CONTAINS
CALL construct_submatrices(matrix_in, subm_in, &
dpattern, map, node_of_domain, select_row)
!!!TRIM IF (matrix_trimmer_required) THEN
!!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
!!!TRIM dpattern,map,node_of_domain,select_row)
!!!TRIM ENDIF
IF (my_action == 0) THEN
! for example, apply preconditioner
CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
@ -2257,16 +2252,10 @@ CONTAINS
ALLOCATE (subm_main(ndomains))
CALL init_submatrices(subm_main)
!!!TRIM ALLOCATE(subm_trimmer(ndomains))
CALL construct_submatrices(matrix_main, subm_main, &
dpattern, map, node_of_domain, select_row_col)
!!!TRIM IF (matrix_trimmer_required) THEN
!!!TRIM CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
!!!TRIM dpattern,map,node_of_domain,select_row)
!!!TRIM ENDIF
IF (my_action == -1) THEN
! project out the local occupied space
!tmp=MATMUL(subm_r(idomain)%mdata,Minv)
@ -2315,38 +2304,9 @@ CONTAINS
END DO
naos = subm_main(idomain)%nrows
!WRITE(*,*) "Domain, mo_self_and_neig, ao_domain: ", idomain, n_domain_mos, naos
ALLOCATE (Minv(naos, naos))
!!!TRIM IF (my_use_trimmer) THEN
!!!TRIM ! THIS IS SUPER EXPENSIVE (ELIMINATE)
!!!TRIM ! trim the main matrix before inverting
!!!TRIM ! assume that the trimmer columns are different (i.e. the main matrix is different for each MO)
!!!TRIM allocate(tmp(naos,nmos(idomain)))
!!!TRIM DO ii=1, nmos(idomain)
!!!TRIM ! transform the main matrix using the trimmer for the current MO
!!!TRIM DO jj=1, naos
!!!TRIM DO kk=1, naos
!!!TRIM Mstore(jj,kk)=sumb_main(idomain)%mdata(jj,kk)*&
!!!TRIM subm_trimmer(idomain)%mdata(jj,ii)*&
!!!TRIM subm_trimmer(idomain)%mdata(kk,ii)
!!!TRIM ENDDO
!!!TRIM ENDDO
!!!TRIM ! invert the main matrix (exclude some eigenvalues, shift some)
!!!TRIM CALL pseudo_invert_matrix(A=Mstore,Ainv=Minv,N=naos,method=1,&
!!!TRIM !range1_thr=1.0E-9_dp,range2_thr=1.0E-9_dp,&
!!!TRIM shift=1.0E-5_dp,&
!!!TRIM range1=nmos(idomain),range2=nmos(idomain),&
!!!TRIM
!!!TRIM ! apply the inverted matrix
!!!TRIM ! RZK-warning this is only possible when the preconditioner is applied
!!!TRIM tmp(:,ii)=MATMUL(Minv,subm_in(idomain)%mdata(:,ii))
!!!TRIM ENDDO
!!!TRIM subm_out=MATMUL(tmp,sigma)
!!!TRIM deallocate(tmp)
!!!TRIM ELSE
IF (PRESENT(bad_modes_projector_down)) THEN
ALLOCATE (proj_array(naos, naos))
CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, &
@ -2360,7 +2320,6 @@ CONTAINS
CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, &
range1=nmos(idomain), range2=n_domain_mos)
END IF
!!!TRIM ENDIF
CALL copy_submatrices(subm_main(idomain), preconditioner(idomain), .FALSE.)
CALL copy_submatrix_data(Minv, preconditioner(idomain))
@ -2635,24 +2594,6 @@ CONTAINS
unit_nr = -1
END IF
!CALL dpotrf('L', N, Ainv, N, INFO )
!IF( INFO/=0 ) THEN
! CPErrorMessage(cp_failure_level,routineP,"DPOTRF failed")
! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
!END IF
!CALL dpotri('L', N, Ainv, N, INFO )
!IF( INFO/=0 ) THEN
! CPErrorMessage(cp_failure_level,routineP,"DPOTRI failed")
! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
!END IF
!! complete the matrix
!DO ii=1,N
! DO jj=ii+1,N
! Ainv(ii,jj)=Ainv(jj,ii)
! ENDDO
! !WRITE(*,'(100F13.9)') Ainv(ii,:)
!ENDDO
! diagonalize first
ALLOCATE (eigenvalues(N))
! Query the optimal workspace for dsyev
@ -2689,21 +2630,6 @@ CONTAINS
DEALLOCATE (eigenvalues)
!!! ! compute the error
!!! allocate(test(N,N))
!!! test=MATMUL(Ainv,A)
!!! DO ii=1,N
!!! test(ii,ii)=test(ii,ii)-1.0_dp
!!! ENDDO
!!! test_error=0.0_dp
!!! DO ii=1,N
!!! DO jj=1,N
!!! test_error=test_error+test(jj,ii)*test(jj,ii)
!!! ENDDO
!!! ENDDO
!!! WRITE(*,*) "Inversion error: ", SQRT(test_error)
!!! deallocate(test)
CALL timestop(handle)
END SUBROUTINE matrix_sqrt
@ -2916,21 +2842,6 @@ CONTAINS
END SELECT
!! compute the inversion error
!allocate(temp1(N,N))
!temp1=MATMUL(Ainv,A)
!DO ii=1,N
! temp1(ii,ii)=temp1(ii,ii)-1.0_dp
!ENDDO
!temp1_error=0.0_dp
!DO ii=1,N
! DO jj=1,N
! temp1_error=temp1_error+temp1(jj,ii)*temp1(jj,ii)
! ENDDO
!ENDDO
!WRITE(*,*) "Inversion error: ", SQRT(temp1_error)
!deallocate(temp1)
CALL timestop(handle)
END SUBROUTINE pseudo_invert_matrix

1640
src/almo_scf_optimizer.F
View file

File diff suppressed because it is too large Load diff

76
src/almo_scf_qs.F
View file

@ -288,63 +288,6 @@ CONTAINS
CALL dbcsr_work_create(matrix_new, work_mutable=.TRUE.)
CALL dbcsr_get_info(matrix_new, nblkrows_total=nblkrows_tot, &
row_blk_size=row_blk_size, col_blk_size=col_blk_size)
! startQQQ - this part of the code scales quadratically
! therefore it is replaced with a less general but linear scaling algorithm below
! the quadratic algorithm is kept to be re-written later
!QQQCALL dbcsr_get_info(matrix_new, nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
!QQQDO row = 1, nblkrows_tot
!QQQ DO col = 1, nblkcols_tot
!QQQ tr = .FALSE.
!QQQ iblock_row = row
!QQQ iblock_col = col
!QQQ CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, tr, hold)
!QQQ IF(hold==mynode) THEN
!QQQ
!QQQ ! RZK-warning replace with a function which says if this
!QQQ ! distribution block is active or not
!QQQ ! Translate indeces of distribution blocks to domain blocks
!QQQ if (size_keys(1)==almo_mat_dim_aobasis) then
!QQQ domain_row=almo_scf_env%domain_index_of_ao_block(iblock_row)
!QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt_full) then
!QQQ domain_row=almo_scf_env%domain_index_of_mo_block(iblock_row)
!QQQ else
!QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension")
!QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
!QQQ endif
!QQQ if (size_keys(2)==almo_mat_dim_aobasis) then
!QQQ domain_col=almo_scf_env%domain_index_of_ao_block(iblock_col)
!QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. &
!QQQ size_keys(2)==almo_mat_dim_virt_full) then
!QQQ domain_col=almo_scf_env%domain_index_of_mo_block(iblock_col)
!QQQ else
!QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension")
!QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
!QQQ endif
!QQQ ! Finds if we need this block
!QQQ ! only the block-diagonal constraint is implemented here
!QQQ active=.false.
!QQQ if (domain_row==domain_col) active=.true.
!QQQ IF (active) THEN
!QQQ ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
!QQQ new_block(:, :) = 1.0_dp
!QQQ CALL dbcsr_put_block(matrix_new, iblock_row, iblock_col, new_block)
!QQQ DEALLOCATE (new_block)
!QQQ ENDIF
!QQQ ENDIF ! mynode
!QQQ ENDDO
!QQQENDDO
!QQQtake care of zero-electron fragments
! endQQQ - end of the quadratic part
! start linear-scaling replacement:
! works only for molecular blocks AND molecular distributions
DO row = 1, nblkrows_tot
@ -611,7 +554,7 @@ CONTAINS
n_el_f=REAL(almo_scf_env%nelectrons_total, dp), &
maxocc=2.0_dp, &
flexible_electron_count=dft_control%relax_multiplicity)
ELSEIF (almo_scf_env%nspins == 2) THEN
ELSE IF (almo_scf_env%nspins == 2) THEN
CALL allocate_mo_set(mo_set=mos(ispin), &
nao=nrow_fm, &
nmo=ncol_fm, &
@ -1496,23 +1439,6 @@ CONTAINS
DEALLOCATE (last_atom_of_molecule)
END IF
!mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(&
! dbcsr_distribution(almo_scf_env%quench_t(ispin))))
!CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
! nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
!DO row = 1, nblkrows_tot
! DO col = 1, nblkcols_tot
! tr = .FALSE.
! iblock_row = row
! iblock_col = col
! CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin),&
! iblock_row, iblock_col, tr, hold)
! CALL dbcsr_get_block_p(almo_scf_env%quench_t(ispin),&
! row, col, p_old_block, found)
! write(*,*) "RST_NOTE:", mynode, row, col, hold, found
! ENDDO
!ENDDO
CALL timestop(handle)
END SUBROUTINE almo_scf_construct_quencher

6
src/almo_scf_types.F
View file

@ -569,8 +569,9 @@ CONTAINS
DO istore = 1, MIN(almo_scf_env%almo_history%istore, almo_scf_env%almo_history%nstore)
CALL dbcsr_release(almo_scf_env%almo_history%matrix_p_up_down(ispin, istore))
END DO
IF (almo_scf_env%almo_history%istore > 0) &
IF (almo_scf_env%almo_history%istore > 0) THEN
CALL dbcsr_release(almo_scf_env%almo_history%matrix_t(ispin))
END IF
END DO
DEALLOCATE (almo_scf_env%almo_history%matrix_p_up_down)
DEALLOCATE (almo_scf_env%almo_history%matrix_t)
@ -580,8 +581,9 @@ CONTAINS
CALL dbcsr_release(almo_scf_env%xalmo_history%matrix_p_up_down(ispin, istore))
!CALL dbcsr_release(almo_scf_env%xalmo_history%matrix_x(ispin, istore))
END DO
IF (almo_scf_env%xalmo_history%istore > 0) &
IF (almo_scf_env%xalmo_history%istore > 0) THEN
CALL dbcsr_release(almo_scf_env%xalmo_history%matrix_t(ispin))
END IF
END DO
DEALLOCATE (almo_scf_env%xalmo_history%matrix_p_up_down)
!DEALLOCATE (almo_scf_env%xalmo_history%matrix_x)

2
src/aobasis/ai_contraction.F
View file

@ -439,7 +439,7 @@ CONTAINS
ELSE
qab(ia:ja, ib:jb) = qab(ia:ja, ib:jb) + sab(1:na, 1:nb)
END IF
ELSEIF (dir == "OUT" .OR. dir == "out") THEN
ELSE IF (dir == "OUT" .OR. dir == "out") THEN
! SAB <= QAB(block)
ja = ia + na - 1
jb = ib + nb - 1

16
src/aobasis/ai_eri_debug.F
View file

@ -128,16 +128,16 @@ CONTAINS
IF (dn(1) > 0) THEN
IABCD = os(an, bn, cn + i1, dn - i1) - (D(1) - C(1))*os(an, bn, cn, dn - i1)
ELSEIF (dn(2) > 0) THEN
ELSE IF (dn(2) > 0) THEN
IABCD = os(an, bn, cn + i2, dn - i2) - (D(2) - C(2))*os(an, bn, cn, dn - i2)
ELSEIF (dn(3) > 0) THEN
ELSE IF (dn(3) > 0) THEN
IABCD = os(an, bn, cn + i3, dn - i3) - (D(3) - C(3))*os(an, bn, cn, dn - i3)
ELSE
IF (bn(1) > 0) THEN
IABCD = os(an + i1, bn - i1, cn, dn) - (B(1) - A(1))*os(an, bn - i1, cn, dn)
ELSEIF (bn(2) > 0) THEN
ELSE IF (bn(2) > 0) THEN
IABCD = os(an + i2, bn - i2, cn, dn) - (B(2) - A(2))*os(an, bn - i2, cn, dn)
ELSEIF (bn(3) > 0) THEN
ELSE IF (bn(3) > 0) THEN
IABCD = os(an + i3, bn - i3, cn, dn) - (B(3) - A(3))*os(an, bn - i3, cn, dn)
ELSE
IF (cn(1) > 0) THEN
@ -145,12 +145,12 @@ CONTAINS
0.5_dp*an(1)/eta*os(an - i1, bn, cn - i1, dn) + &
0.5_dp*(cn(1) - 1)/eta*os(an, bn, cn - i1 - i1, dn) - &
xsi/eta*os(an + i1, bn, cn - i1, dn)
ELSEIF (cn(2) > 0) THEN
ELSE IF (cn(2) > 0) THEN
IABCD = ((Q(2) - C(2)) + xsi/eta*(P(2) - A(2)))*os(an, bn, cn - i2, dn) + &
0.5_dp*an(2)/eta*os(an - i2, bn, cn - i2, dn) + &
0.5_dp*(cn(2) - 1)/eta*os(an, bn, cn - i2 - i2, dn) - &
xsi/eta*os(an + i2, bn, cn - i2, dn)
ELSEIF (cn(3) > 0) THEN
ELSE IF (cn(3) > 0) THEN
IABCD = ((Q(3) - C(3)) + xsi/eta*(P(3) - A(3)))*os(an, bn, cn - i3, dn) + &
0.5_dp*an(3)/eta*os(an - i3, bn, cn - i3, dn) + &
0.5_dp*(cn(3) - 1)/eta*os(an, bn, cn - i3 - i3, dn) - &
@ -161,12 +161,12 @@ CONTAINS
(W(1) - P(1))*os(an - i1, bn, cn, dn, m + 1) + &
0.5_dp*(an(1) - 1)/xsi*os(an - i1 - i1, bn, cn, dn, m) - &
0.5_dp*(an(1) - 1)/xsi*rho/xsi*os(an - i1 - i1, bn, cn, dn, m + 1)
ELSEIF (an(2) > 0) THEN
ELSE IF (an(2) > 0) THEN
IABCD = (P(2) - A(2))*os(an - i2, bn, cn, dn, m) + &
(W(2) - P(2))*os(an - i2, bn, cn, dn, m + 1) + &
0.5_dp*(an(2) - 1)/xsi*os(an - i2 - i2, bn, cn, dn, m) - &
0.5_dp*(an(2) - 1)/xsi*rho/xsi*os(an - i2 - i2, bn, cn, dn, m + 1)
ELSEIF (an(3) > 0) THEN
ELSE IF (an(3) > 0) THEN
IABCD = (P(3) - A(3))*os(an - i3, bn, cn, dn, m) + &
(W(3) - P(3))*os(an - i3, bn, cn, dn, m + 1) + &
0.5_dp*(an(3) - 1)/xsi*os(an - i3 - i3, bn, cn, dn, m) - &

8
src/aobasis/ai_os_rr.F
View file

@ -162,7 +162,7 @@ CONTAINS
rr(m, coa, 1) = rr(m, coa, 1) + g*REAL(az - 1, dp)*(rr(m, coa2z, 1) - rr(m + 1, coa2z, 1))
END DO
END IF
ELSEIF (ay > 0) THEN
ELSE IF (ay > 0) THEN
DO m = 0, mmax - la
rr(m, coa, 1) = rap(2)*rr(m, coa1y, 1) - rcp(2)*rr(m + 1, coa1y, 1)
END DO
@ -171,7 +171,7 @@ CONTAINS
rr(m, coa, 1) = rr(m, coa, 1) + g*REAL(ay - 1, dp)*(rr(m, coa2y, 1) - rr(m + 1, coa2y, 1))
END DO
END IF
ELSEIF (ax > 0) THEN
ELSE IF (ax > 0) THEN
DO m = 0, mmax - la
rr(m, coa, 1) = rap(1)*rr(m, coa1x, 1) - rcp(1)*rr(m + 1, coa1x, 1)
END DO
@ -225,7 +225,7 @@ CONTAINS
rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(az, dp)*(rr(m, coa1z, cob1z) - rr(m + 1, coa1z, cob1z))
END DO
END IF
ELSEIF (by > 0) THEN
ELSE IF (by > 0) THEN
DO m = 0, mmax - la - lb
rr(m, coa, cob) = rbp(2)*rr(m, coa, cob1y) - rcp(2)*rr(m + 1, coa, cob1y)
END DO
@ -239,7 +239,7 @@ CONTAINS
rr(m, coa, cob) = rr(m, coa, cob) + g*REAL(ay, dp)*(rr(m, coa1y, cob1y) - rr(m + 1, coa1y, cob1y))
END DO
END IF
ELSEIF (bx > 0) THEN
ELSE IF (bx > 0) THEN
DO m = 0, mmax - la - lb
rr(m, coa, cob) = rbp(1)*rr(m, coa, cob1x) - rcp(1)*rr(m + 1, coa, cob1x)
END DO

12
src/aobasis/ai_overlap.F
View file

@ -2050,12 +2050,12 @@ CONTAINS
DO l = 0, la
IF (l == 0) THEN
fun(:, l) = z_one
ELSEIF (l == 1) THEN
ELSE IF (l == 1) THEN
fun(:, l) = CMPLX(0.0_dp, 0.5_dp*oa*gval(:), KIND=dp)
ELSEIF (l == 2) THEN
ELSE IF (l == 2) THEN
fun(:, l) = CMPLX(-(0.5_dp*oa*gval(:))**2, 0.0_dp, KIND=dp)
fun(:, l) = fun(:, l) + CMPLX(0.5_dp*oa, 0.0_dp, KIND=dp)
ELSEIF (l == 3) THEN
ELSE IF (l == 3) THEN
fun(:, l) = CMPLX(0.0_dp, -(0.5_dp*oa*gval(:))**3, KIND=dp)
fun(:, l) = fun(:, l) + CMPLX(0.0_dp, 0.75_dp*oa*oa*gval(:), KIND=dp)
ELSE
@ -2065,12 +2065,12 @@ CONTAINS
DO l = 0, lb
IF (l == 0) THEN
gun(:, l) = z_one
ELSEIF (l == 1) THEN
ELSE IF (l == 1) THEN
gun(:, l) = CMPLX(0.0_dp, 0.5_dp*ob*gval(:), KIND=dp)
ELSEIF (l == 2) THEN
ELSE IF (l == 2) THEN
gun(:, l) = CMPLX(-(0.5_dp*ob*gval(:))**2, 0.0_dp, KIND=dp)
gun(:, l) = gun(:, l) + CMPLX(0.5_dp*ob, 0.0_dp, KIND=dp)
ELSEIF (l == 3) THEN
ELSE IF (l == 3) THEN
gun(:, l) = CMPLX(0.0_dp, -(0.5_dp*ob*gval(:))**3, KIND=dp)
gun(:, l) = gun(:, l) + CMPLX(0.0_dp, 0.75_dp*ob*ob*gval(:), KIND=dp)
ELSE

12
src/aobasis/ai_overlap3_debug.F
View file

@ -99,25 +99,25 @@ CONTAINS
IF (bn(1) > 0) THEN
IACB = os_overlap3(an, cn + i1, bn - i1) + (C(1) - B(1))*os_overlap3(an, cn, bn - i1)
ELSEIF (bn(2) > 0) THEN
ELSE IF (bn(2) > 0) THEN
IACB = os_overlap3(an, cn + i2, bn - i2) + (C(2) - B(2))*os_overlap3(an, cn, bn - i2)
ELSEIF (bn(3) > 0) THEN
ELSE IF (bn(3) > 0) THEN
IACB = os_overlap3(an, cn + i3, bn - i3) + (C(3) - B(3))*os_overlap3(an, cn, bn - i3)
ELSE
IF (cn(1) > 0) THEN
IACB = os_overlap3(an + i1, cn - i1, bn) + (A(1) - C(1))*os_overlap3(an, cn - i1, bn)
ELSEIF (cn(2) > 0) THEN
ELSE IF (cn(2) > 0) THEN
IACB = os_overlap3(an + i2, cn - i2, bn) + (A(2) - C(2))*os_overlap3(an, cn - i2, bn)
ELSEIF (cn(3) > 0) THEN
ELSE IF (cn(3) > 0) THEN
IACB = os_overlap3(an + i3, cn - i3, bn) + (A(3) - C(3))*os_overlap3(an, cn - i3, bn)
ELSE
IF (an(1) > 0) THEN
IACB = (G(1) - A(1))*os_overlap3(an - i1, cn, bn) + &
0.5_dp*(an(1) - 1)/(xsi + xc)*os_overlap3(an - i1 - i1, cn, bn)
ELSEIF (an(2) > 0) THEN
ELSE IF (an(2) > 0) THEN
IACB = (G(2) - A(2))*os_overlap3(an - i2, cn, bn) + &
0.5_dp*(an(2) - 1)/(xsi + xc)*os_overlap3(an - i2 - i2, cn, bn)
ELSEIF (an(3) > 0) THEN
ELSE IF (an(3) > 0) THEN
IACB = (G(3) - A(3))*os_overlap3(an - i3, cn, bn) + &
0.5_dp*(an(3) - 1)/(xsi + xc)*os_overlap3(an - i3 - i3, cn, bn)
ELSE

8
src/aobasis/ai_overlap_debug.F
View file

@ -87,18 +87,18 @@ CONTAINS
IF (bn(1) > 0) THEN
IAB = os_overlap2(an + i1, bn - i1) + (A(1) - B(1))*os_overlap2(an, bn - i1)
ELSEIF (bn(2) > 0) THEN
ELSE IF (bn(2) > 0) THEN
IAB = os_overlap2(an + i2, bn - i2) + (A(2) - B(2))*os_overlap2(an, bn - i2)
ELSEIF (bn(3) > 0) THEN
ELSE IF (bn(3) > 0) THEN
IAB = os_overlap2(an + i3, bn - i3) + (A(3) - B(3))*os_overlap2(an, bn - i3)
ELSE
IF (an(1) > 0) THEN
IAB = (P(1) - A(1))*os_overlap2(an - i1, bn) + &
0.5_dp*(an(1) - 1)/xsi*os_overlap2(an - i1 - i1, bn)
ELSEIF (an(2) > 0) THEN
ELSE IF (an(2) > 0) THEN
IAB = (P(2) - A(2))*os_overlap2(an - i2, bn) + &
0.5_dp*(an(2) - 1)/xsi*os_overlap2(an - i2 - i2, bn)
ELSEIF (an(3) > 0) THEN
ELSE IF (an(3) > 0) THEN
IAB = (P(3) - A(3))*os_overlap2(an - i3, bn) + &
0.5_dp*(an(3) - 1)/xsi*os_overlap2(an - i3 - i3, bn)
ELSE

8
src/aobasis/basis_set_types.F
View file

@ -1680,8 +1680,9 @@ CONTAINS
l(nshell(iset) - ishell + i, iset) = lshell
END DO
END DO
IF (LEN_TRIM(line_att) /= 0) &
IF (LEN_TRIM(line_att) /= 0) THEN
CPABORT("Error reading the Basis from input file!")
END IF
DO ipgf = 1, npgf(iset)
is_ok = cp_sll_val_next(list, val)
IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!")
@ -2503,8 +2504,9 @@ CONTAINS
ng = gto_basis_set%npgf(1)
DO iset = 1, nset
IF ((ng /= gto_basis_set%npgf(iset)) .AND. do_ortho) &
IF ((ng /= gto_basis_set%npgf(iset)) .AND. do_ortho) THEN
CPABORT("different number of primitves")
END IF
END DO
IF (do_ortho) THEN
@ -2682,7 +2684,7 @@ CONTAINS
s00 = ai*aj*(pi*ab)**1.50_dp
IF (l == 0) THEN
sss = s00
ELSEIF (l == 1) THEN
ELSE IF (l == 1) THEN
sss = s00*ab*0.5_dp
ELSE
CPABORT("aovlp lvalue")

27
src/arnoldi/arnoldi_data_methods.F
View file

@ -138,13 +138,15 @@ CONTAINS
CALL control%mp_group%set_handle(group_handle)
CALL control%pcol_group%set_handle(pcol_handle)
IF (.NOT. subgroups_defined) &
IF (.NOT. subgroups_defined) THEN
CPABORT("arnoldi only with subgroups")
END IF
control%symmetric = .FALSE.
! Will need a fix for complex because there it has to be hermitian
IF (SIZE(matrix) == 1) &
IF (SIZE(matrix) == 1) THEN
control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
END IF
! Set the control parameters
control%max_iter = max_iter
@ -158,23 +160,28 @@ CONTAINS
control%nrestart = nrestarts
control%generalized_ev = generalized_ev
IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) &
IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) THEN
CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// &
'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0')
END IF
! some checks for the generalized EV mode
IF (control%generalized_ev .AND. selection_crit == 1) &
IF (control%generalized_ev .AND. selection_crit == 1) THEN
CALL cp_abort(__LOCATION__, &
'generalized ev can only highest OR lowest EV')
IF (control%generalized_ev .AND. nval_request /= 1) &
END IF
IF (control%generalized_ev .AND. nval_request /= 1) THEN
CALL cp_abort(__LOCATION__, &
'generalized ev can only compute one EV at the time')
IF (control%generalized_ev .AND. control%nrestart == 0) &
END IF
IF (control%generalized_ev .AND. control%nrestart == 0) THEN
CALL cp_abort(__LOCATION__, &
'outer loops are mandatory for generalized EV, set nrestart appropriatly')
IF (SIZE(matrix) /= 2 .AND. control%generalized_ev) &
END IF
IF (SIZE(matrix) /= 2 .AND. control%generalized_ev) THEN
CALL cp_abort(__LOCATION__, &
'generalized ev needs exactly two matrices as input (2nd is the metric)')
END IF
ALLOCATE (control%selected_ind(max_iter))
CALL set_control(arnoldi_env, control)
@ -386,8 +393,9 @@ CONTAINS
INTEGER :: ev_ind
INTEGER, DIMENSION(:), POINTER :: selected_ind
IF (ind > get_nval_out(arnoldi_env)) &
IF (ind > get_nval_out(arnoldi_env)) THEN
CPABORT('outside range of indexed evals')
END IF
selected_ind => get_sel_ind(arnoldi_env)
ev_ind = selected_ind(ind)
@ -411,8 +419,9 @@ CONTAINS
INTEGER, DIMENSION(:), POINTER :: selected_ind
NULLIFY (evals)
IF (SIZE(eval_out) < get_nval_out(arnoldi_env)) &
IF (SIZE(eval_out) < get_nval_out(arnoldi_env)) THEN
CPABORT('array for eval output too small')
END IF
selected_ind => get_sel_ind(arnoldi_env)
evals => get_evals(arnoldi_env)

2
src/arnoldi/arnoldi_geev.F
View file

@ -144,7 +144,7 @@ CONTAINS
i = 1
DO WHILE (i <= ndim)
IF (ABS(eval2(i)) < EPSILON(REAL(0.0, dp))) THEN
evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i)))
evec_r(:, i) = evec_r(:, i)/NORM2(evec_r(:, i))
revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, dp), dp)
i = i + 1

2
src/arnoldi/arnoldi_methods.F
View file

@ -561,7 +561,7 @@ CONTAINS
ar_data%local_history = Zmat
! broadcast the Hessenberg matrix so we don't need to care later on
DEALLOCATE (v_vec); DEALLOCATE (w_vec); DEALLOCATE (s_vec); DEALLOCATE (h_vec); DEALLOCATE (CZmat);
DEALLOCATE (v_vec); DEALLOCATE (w_vec); DEALLOCATE (s_vec); DEALLOCATE (h_vec); DEALLOCATE (CZmat)
DEALLOCATE (Zmat); DEALLOCATE (BZmat)
CALL timestop(handle)

8
src/atom_electronic_structure.F
View file

@ -343,10 +343,10 @@ CONTAINS
xcmat%op = 0._dp
CALL calculate_atom_vxc_lda(xcmat, atom, xc_section)
! ZMP added options for the zmp calculations, building external density and vxc potential
ELSEIF (need_zmp) THEN
ELSE IF (need_zmp) THEN
xcmat%op = 0._dp
CALL calculate_atom_zmp(ext_density=ext_density, atom=atom, lprint=.FALSE., xcmat=xcmat)
ELSEIF (need_vxc) THEN
ELSE IF (need_vxc) THEN
xcmat%op = 0._dp
CALL calculate_atom_ext_vxc(vxc=ext_vxc, atom=atom, lprint=.FALSE., xcmat=xcmat)
ELSE
@ -503,10 +503,10 @@ CONTAINS
ne = atom%state%occupation(l, k)
IF (ne == 0._dp) THEN !empty shell
EXIT !assume there are no holes
ELSEIF (ne == 2._dp*nm) THEN !closed shell
ELSE IF (ne == 2._dp*nm) THEN !closed shell
atom%state%occa(l, k) = nm
atom%state%occb(l, k) = nm
ELSEIF (atom%state%multiplicity == -2) THEN !High spin case
ELSE IF (atom%state%multiplicity == -2) THEN !High spin case
atom%state%occa(l, k) = MIN(ne, nm)
atom%state%occb(l, k) = MAX(0._dp, ne - nm)
ELSE

6
src/atom_energy.F
View file

@ -1103,11 +1103,11 @@ CONTAINS
IF (PRESENT(counter)) THEN
WRITE (str, "(I12)") counter
ELSEIF (PRESENT(rval)) THEN
ELSE IF (PRESENT(rval)) THEN
WRITE (str, "(G18.8)") rval
ELSEIF (PRESENT(ival)) THEN
ELSE IF (PRESENT(ival)) THEN
WRITE (str, "(I12)") ival
ELSEIF (PRESENT(cval)) THEN
ELSE IF (PRESENT(cval)) THEN
WRITE (str, "(A)") TRIM(ADJUSTL(cval))
ELSE
WRITE (str, "(A)") ""

4
src/atom_fit.F
View file

@ -657,7 +657,7 @@ CONTAINS
ntarget = ntarget + 1
wtot = wtot + atom%weight*w_virt/100._dp
END IF
ELSEIF (k < atom%state%maxn_occ(l)) THEN
ELSE IF (k < atom%state%maxn_occ(l)) THEN
atom%orbitals%wrefene(k, l, 1) = w_semi
atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1) = t_semi
@ -735,7 +735,7 @@ CONTAINS
wtot = wtot + atom%weight*2._dp*w_virt/100._dp
ntarget = ntarget + 2
END IF
ELSEIF (k < atom%state%maxn_occ(l)) THEN
ELSE IF (k < atom%state%maxn_occ(l)) THEN
atom%orbitals%wrefene(k, l, 1:2) = w_semi
atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp
atom%orbitals%crefene(k, l, 1:2) = t_semi

3
src/atom_grb.F
View file

@ -152,8 +152,9 @@ CONTAINS
CALL allocate_grid_atom(basis%grid)
CALL section_vals_val_get(grb_section, "QUADRATURE", i_val=quadtype)
CALL section_vals_val_get(grb_section, "GRID_POINTS", i_val=ngp)
IF (ngp <= 0) &
IF (ngp <= 0) THEN
CPABORT("# point radial grid < 0")
END IF
CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
basis%grid%nr = ngp
!

2
src/atom_kind_orbitals.F
View file

@ -273,7 +273,7 @@ CONTAINS
! total number of occupied orbitals
IF (PRESENT(nocc) .AND. ghost) THEN
nocc = 0
ELSEIF (PRESENT(nocc)) THEN
ELSE IF (PRESENT(nocc)) THEN
nocc = 0
DO l = 0, lmat
DO k = 1, 7

4
src/atom_operators.F
View file

@ -120,14 +120,14 @@ CONTAINS
rc = potential%rcon
sc = potential%scon
cpot(1:m) = (basis%grid%rad(1:m)/rc)**sc
ELSEIF (potential%conf_type == barrier_conf) THEN
ELSE IF (potential%conf_type == barrier_conf) THEN
om = potential%rcon
ron = potential%scon
rc = ron + om
DO i = 1, m
IF (basis%grid%rad(i) < ron) THEN
cpot(i) = 0.0_dp
ELSEIF (basis%grid%rad(i) < rc) THEN
ELSE IF (basis%grid%rad(i) < rc) THEN
x = (basis%grid%rad(i) - ron)/om
x = 1._dp - x
cpot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp

2
src/atom_optimization.F
View file

@ -221,7 +221,7 @@ CONTAINS
IF (nm < 1) nm = history%max_history
fmat = a*history%hmat(nnow)%fmat + (1._dp - a)*history%hmat(nm)%fmat
END IF
ELSEIF (history%hlen == 1) THEN
ELSE IF (history%hlen == 1) THEN
fmat = history%hmat(nnow)%fmat
ELSE
CPABORT("Length of matrix history hlen < 1")

9
src/atom_output.F
View file

@ -192,16 +192,19 @@ CONTAINS
WRITE (iw, '(T36,A,T61,F20.12)') " Virial (-V/T) ::", -atom%energy%epot/atom%energy%ekin
END IF
WRITE (iw, '(T36,A,T61,F20.12)') " Core Energy ::", atom%energy%ecore
IF (atom%energy%exc /= 0._dp) &
IF (atom%energy%exc /= 0._dp) THEN
WRITE (iw, '(T36,A,T61,F20.12)') " XC Energy ::", atom%energy%exc
END IF
WRITE (iw, '(T36,A,T61,F20.12)') " Coulomb Energy ::", atom%energy%ecoulomb
IF (atom%energy%eexchange /= 0._dp) &
IF (atom%energy%eexchange /= 0._dp) THEN
WRITE (iw, '(T34,A,T61,F20.12)') "HF Exchange Energy ::", atom%energy%eexchange
END IF
IF (atom%potential%ppot_type /= NO_PSEUDO) THEN
WRITE (iw, '(T20,A,T61,F20.12)') " Total Pseudopotential Energy ::", atom%energy%epseudo
WRITE (iw, '(T20,A,T61,F20.12)') " Local Pseudopotential Energy ::", atom%energy%eploc
IF (atom%energy%elsd /= 0._dp) &
IF (atom%energy%elsd /= 0._dp) THEN
WRITE (iw, '(T20,A,T61,F20.12)') " Local Spin-potential Energy ::", atom%energy%elsd
END IF
WRITE (iw, '(T20,A,T61,F20.12)') " Nonlocal Pseudopotential Energy ::", atom%energy%epnl
END IF
IF (atom%potential%confinement) THEN

4
src/atom_pseudo.F
View file

@ -257,10 +257,10 @@ CONTAINS
ne = state%occupation(l, k)
IF (ne == 0._dp) THEN !empty shell
EXIT !assume there are no holes
ELSEIF (ne == 2._dp*nm) THEN !closed shell
ELSE IF (ne == 2._dp*nm) THEN !closed shell
state%occa(l, k) = nm
state%occb(l, k) = nm
ELSEIF (state%multiplicity == -2) THEN !High spin case
ELSE IF (state%multiplicity == -2) THEN !High spin case
state%occa(l, k) = MIN(ne, nm)
state%occb(l, k) = MAX(0._dp, ne - nm)
ELSE

2
src/atom_set_basis.F
View file

@ -204,7 +204,7 @@ CONTAINS
basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
END DO
ELSEIF (basis%basis_type == CGTO_BASIS) THEN
ELSE IF (basis%basis_type == CGTO_BASIS) THEN
DO k = 1, nr
rk = basis%grid%rad(k)
ear = EXP(-al*basis%grid%rad(k)**2)

8
src/atom_sgp.F
View file

@ -122,7 +122,7 @@ CONTAINS
! generate the transformed potentials
IF (is_ecp) THEN
CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
ELSEIF (is_upf) THEN
ELSE IF (is_upf) THEN
CALL upf_sgp_constr(upf_pot, sgp_pot, basis)
ELSE
CPABORT("Either ecp_pot or upf_pot is needed for sgp_construction")
@ -137,7 +137,7 @@ CONTAINS
!
IF (is_ecp) THEN
CALL ecpints(hnl%op, basis, ecp_pot)
ELSEIF (is_upf) THEN
ELSE IF (is_upf) THEN
CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
ELSE
CPABORT("Either ecp_pot or upf_pot is needed for sgp_construction")
@ -246,7 +246,7 @@ CONTAINS
IF (do_transform) THEN
IF (is_ecp) THEN
CALL ecp_sgp_constr(ecp_pot, sgp_pot, atom_ref%basis)
ELSEIF (is_upf) THEN
ELSE IF (is_upf) THEN
CALL upf_sgp_constr(upf_pot, sgp_pot, atom_ref%basis)
ELSE
CPABORT("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")
@ -279,7 +279,7 @@ CONTAINS
!
IF (is_ecp) THEN
CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
ELSEIF (is_upf) THEN
ELSE IF (is_upf) THEN
CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
ELSE
CPABORT("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")

6
src/atom_types.F
View file

@ -414,8 +414,9 @@ CONTAINS
CALL allocate_grid_atom(basis%grid)
CALL section_vals_val_get(basis_section, "QUADRATURE", i_val=quadtype)
CALL section_vals_val_get(basis_section, "GRID_POINTS", i_val=ngp)
IF (ngp <= 0) &
IF (ngp <= 0) THEN
CPABORT("The number of radial grid points must be greater than zero.")
END IF
CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
basis%grid%nr = ngp
basis%geometrical = .FALSE.
@ -840,8 +841,9 @@ CONTAINS
CALL allocate_grid_atom(gbasis%grid)
ngp = SIZE(r)
quadtype = do_gapw_log
IF (ngp <= 0) &
IF (ngp <= 0) THEN
CPABORT("The number of radial grid points must be greater than zero.")
END IF
CALL create_grid_atom(gbasis%grid, ngp, 1, 1, 0, quadtype)
gbasis%grid%nr = ngp
gbasis%grid%rad(:) = r(:)

22
src/atom_upf.F
View file

@ -218,39 +218,39 @@ CONTAINS
IF (nametag(2:8) == "PP_INFO") THEN
CPASSERT(nametag(9:9) == ">")
CALL upf_info_section(parser, pot)
ELSEIF (nametag(2:10) == "PP_HEADER") THEN
ELSE IF (nametag(2:10) == "PP_HEADER") THEN
IF (.NOT. (nametag(11:11) == ">")) THEN
CALL upf_header_option(parser, pot)
END IF
ELSEIF (nametag(2:8) == "PP_MESH") THEN
ELSE IF (nametag(2:8) == "PP_MESH") THEN
IF (.NOT. (nametag(9:9) == ">")) THEN
CALL upf_mesh_option(parser, pot)
END IF
CALL upf_mesh_section(parser, pot)
ELSEIF (nametag(2:8) == "PP_NLCC") THEN
ELSE IF (nametag(2:8) == "PP_NLCC") THEN
IF (nametag(9:9) == ">") THEN
CALL upf_nlcc_section(parser, pot, .FALSE.)
ELSE
CALL upf_nlcc_section(parser, pot, .TRUE.)
END IF
ELSEIF (nametag(2:9) == "PP_LOCAL") THEN
ELSE IF (nametag(2:9) == "PP_LOCAL") THEN
IF (nametag(10:10) == ">") THEN
CALL upf_local_section(parser, pot, .FALSE.)
ELSE
CALL upf_local_section(parser, pot, .TRUE.)
END IF
ELSEIF (nametag(2:12) == "PP_NONLOCAL") THEN
ELSE IF (nametag(2:12) == "PP_NONLOCAL") THEN
CPASSERT(nametag(13:13) == ">")
CALL upf_nonlocal_section(parser, pot)
ELSEIF (nametag(2:13) == "PP_SEMILOCAL") THEN
ELSE IF (nametag(2:13) == "PP_SEMILOCAL") THEN
CALL upf_semilocal_section(parser, pot)
ELSEIF (nametag(2:9) == "PP_PSWFC") THEN
ELSE IF (nametag(2:9) == "PP_PSWFC") THEN
! skip section for now
ELSEIF (nametag(2:11) == "PP_RHOATOM") THEN
ELSE IF (nametag(2:11) == "PP_RHOATOM") THEN
! skip section for now
ELSEIF (nametag(2:7) == "PP_PAW") THEN
ELSE IF (nametag(2:7) == "PP_PAW") THEN
! skip section for now
ELSEIF (nametag(2:6) == "/UPF>") THEN
ELSE IF (nametag(2:6) == "/UPF>") THEN
EXIT
END IF
END IF
@ -856,7 +856,7 @@ CONTAINS
END IF
IF (icount > ms) EXIT
END DO
ELSEIF (string(1:15) == "</PP_SEMILOCAL>") THEN
ELSE IF (string(1:15) == "</PP_SEMILOCAL>") THEN
EXIT
ELSE
!

4
src/atom_utils.F
View file

@ -2549,7 +2549,7 @@ CONTAINS
ja = ibptr(ia, la)
jb = ibptr(ib, lb)
smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + sab(1:nna, 1:nnb)
ELSEIF (basis%basis_type == CGTO_BASIS) THEN
ELSE IF (basis%basis_type == CGTO_BASIS) THEN
DO ka = 1, basis%nbas(la)
DO kb = 1, basis%nbas(lb)
ja = ibptr(ka, la)
@ -2572,7 +2572,7 @@ CONTAINS
jb = ibptr(ib, lb)
smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) &
+ sab(1:nna, 1:nnb)
ELSEIF (basis%basis_type == CGTO_BASIS) THEN
ELSE IF (basis%basis_type == CGTO_BASIS) THEN
DO ka = 1, basis%nbas(la)
DO kb = 1, basis%nbas(lb)
ja = ibptr(ka, la)

15
src/atoms_input.F
View file

@ -144,8 +144,9 @@ CONTAINS
EXIT
END IF
END DO
IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) &
IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) THEN
CPABORT("incorrectly formatted line in coord section'"//line_att//"'")
END IF
IF (wrd == 1) THEN
atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1)))
ELSE
@ -195,24 +196,26 @@ CONTAINS
EXIT
END IF
END DO
IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) &
IF (LEN_TRIM(line_att(start_c:end_c - 1)) == 0) THEN
CALL cp_abort(__LOCATION__, &
"Incorrectly formatted input line for atom "// &
TRIM(ADJUSTL(cp_to_string(iatom)))// &
" found in COORD section. Input line: <"// &
TRIM(line_att)//"> ")
END IF
SELECT CASE (wrd)
CASE (1)
atom_info%id_atmname(iatom) = str2id(s2s(line_att(start_c:end_c - 1)))
CASE (2:4)
CALL read_float_object(line_att(start_c:end_c - 1), &
atom_info%r(wrd - 1, iatom), error_message)
IF (LEN_TRIM(error_message) /= 0) &
IF (LEN_TRIM(error_message) /= 0) THEN
CALL cp_abort(__LOCATION__, &
"Incorrectly formatted input line for atom "// &
TRIM(ADJUSTL(cp_to_string(iatom)))// &
" found in COORD section. "//TRIM(error_message)// &
" Input line: <"//TRIM(line_att)//"> ")
END IF
CASE (5)
READ (line_att(start_c:end_c - 1), *) strtmp
atom_info%id_molname(iatom) = str2id(strtmp)
@ -343,8 +346,9 @@ CONTAINS
EXIT
END IF
END DO
IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) &
IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) THEN
CPABORT("incorrectly formatted line in coord section'"//line_att//"'")
END IF
IF (wrd == 1) THEN
at_name(ishell) = line_att(start_c:end_c - 1)
CALL uppercase(at_name(ishell))
@ -393,8 +397,9 @@ CONTAINS
EXIT
END IF
END DO
IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) &
IF (wrd /= 5 .AND. end_c >= LEN(line_att) + 1) THEN
CPABORT("incorrectly formatted line in coord section'"//line_att//"'")
END IF
IF (wrd == 1) THEN
at_name_c(ishell) = line_att(start_c:end_c - 1)
CALL uppercase(at_name_c(ishell))

3
src/base/base_hooks.F
View file

@ -145,8 +145,9 @@ CONTAINS
IF (ASSOCIATED(timestop_hook)) THEN
CALL timestop_hook(handle)
ELSE
IF (handle /= -1) &
IF (handle /= -1) THEN
CALL cp_abort(cp__l("base_hooks.F", __LINE__), "Got wrong handle")
END IF
END IF
END SUBROUTINE timestop

9
src/base/machine.F
View file

@ -741,14 +741,17 @@ CONTAINS
! on a posix system LOGNAME should be defined
CALL get_environment_variable("LOGNAME", value=user, status=istat)
! nope, check alternative
IF (istat /= 0) &
IF (istat /= 0) THEN
CALL get_environment_variable("USER", value=user, status=istat)
END IF
! nope, check alternative
IF (istat /= 0) &
IF (istat /= 0) THEN
CALL get_environment_variable("USERNAME", value=user, status=istat)
END IF
! fall back
IF (istat /= 0) &
IF (istat /= 0) THEN
user = "<unknown>"
END IF
END SUBROUTINE m_getlog

4
src/bse_util.F
View file

@ -1481,7 +1481,7 @@ CONTAINS
cell, dft_control, particle_set, pw_env)
IF (iset == 1) THEN
WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State"
ELSEIF (iset == 2) THEN
ELSE IF (iset == 2) THEN
WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State"
END IF
info_approx_trunc = TRIM(ADJUSTL(info_approximation))
@ -1493,7 +1493,7 @@ CONTAINS
log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io)
IF (iset == 1) THEN
WRITE (title, *) "Natural Transition Orbital Hole State", i
ELSEIF (iset == 2) THEN
ELSE IF (iset == 2) THEN
WRITE (title, *) "Natural Transition Orbital Particle State", i
END IF
CALL cp_pw_to_cube(wf_r, unit_nr_cube, title, particles=particles, stride=stride, mpi_io=mpi_io)

10
src/bsse.F
View file

@ -370,10 +370,10 @@ CONTAINS
iw = cp_print_key_unit_nr(logger, bsse_section, "PRINT%PROGRAM_RUN_INFO", &
extension=".log")
IF (iw > 0) THEN
WRITE (conf_s, fmt="(1000I0)", iostat=istat) conf;
WRITE (conf_s, fmt="(1000I0)", iostat=istat) conf
IF (istat /= 0) conf_s = "exceeded"
CALL compress(conf_s, full=.TRUE.)
WRITE (conf_loc_s, fmt="(1000I0)", iostat=istat) conf_loc;
WRITE (conf_loc_s, fmt="(1000I0)", iostat=istat) conf_loc
IF (istat /= 0) conf_loc_s = "exceeded"
CALL compress(conf_loc_s, full=.TRUE.)
@ -429,15 +429,17 @@ CONTAINS
IF (explicit) THEN
DO i = 1, nconf
CALL section_vals_val_get(configurations, "GLB_CONF", i_rep_section=i, i_vals=glb_conf)
IF (SIZE(glb_conf) /= SIZE(conf)) &
IF (SIZE(glb_conf) /= SIZE(conf)) THEN
CALL cp_abort(__LOCATION__, &
"GLB_CONF requires a binary description of the configuration. Number of integer "// &
"different from the number of fragments defined!")
END IF
CALL section_vals_val_get(configurations, "SUB_CONF", i_rep_section=i, i_vals=sub_conf)
IF (SIZE(sub_conf) /= SIZE(conf)) &
IF (SIZE(sub_conf) /= SIZE(conf)) THEN
CALL cp_abort(__LOCATION__, &
"SUB_CONF requires a binary description of the configuration. Number of integer "// &
"different from the number of fragments defined!")
END IF
IF (ALL(conf == glb_conf) .AND. ALL(conf_loc == sub_conf)) THEN
CALL section_vals_val_get(configurations, "CHARGE", i_rep_section=i, &
i_val=present_charge)

33
src/cell_methods.F
View file

@ -461,11 +461,12 @@ CONTAINS
CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
IF (cell_read_file) THEN ! Case 1
tmp_comb_cell = (cell_read_abc .OR. (cell_read_a .OR. (cell_read_b .OR. cell_read_c)))
IF (tmp_comb_cell) &
IF (tmp_comb_cell) THEN
CALL cp_warn(__LOCATION__, &
"Cell Information provided through A, B, C, or ABC in conjunction "// &
"with CELL_FILE_NAME. The definition in external file will override "// &
"other ones.")
END IF
CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=cell_file_format)
SELECT CASE (cell_file_format)
@ -491,10 +492,11 @@ CONTAINS
read_len = cell_par
CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_par)
read_ang = cell_par
IF (cell_read_a .OR. cell_read_b .OR. cell_read_c) &
IF (cell_read_a .OR. cell_read_b .OR. cell_read_c) THEN
CALL cp_warn(__LOCATION__, &
"Cell information provided through vectors A, B or C in conjunction with ABC. "// &
"The definition of the ABC keyword will override the one provided by A, B and C.")
END IF
ELSE ! Case 3
tmp_comb_abc = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
IF (tmp_comb_abc) THEN
@ -504,10 +506,11 @@ CONTAINS
read_mat(:, 2) = cell_par(:)
CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
read_mat(:, 3) = cell_par(:)
IF (cell_read_alpha_beta_gamma) &
IF (cell_read_alpha_beta_gamma) THEN
CALL cp_warn(__LOCATION__, &
"The keyword ALPHA_BETA_GAMMA is ignored because it was used without the "// &
"keyword ABC.")
END IF
ELSE
CALL cp_abort(__LOCATION__, &
"Neither of the keywords CELL_FILE_NAME or ABC are specified, "// &
@ -722,10 +725,11 @@ CONTAINS
CPASSERT(ASSOCIATED(cell))
! Abort, if one of the value is set to zero
IF (ANY(multiple_unit_cell <= 0)) &
IF (ANY(multiple_unit_cell <= 0)) THEN
CALL cp_abort(__LOCATION__, &
"CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
"A value of 0 or negative is meaningless!")
END IF
! Scale abc according to user request
cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
@ -778,8 +782,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_a", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_lengths(1))
cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
@ -790,8 +795,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_b", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_lengths(2))
cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
@ -802,8 +808,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.length_c", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_lengths(3))
cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
@ -814,8 +821,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_angles(1))
cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
@ -826,8 +834,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_angles(2))
cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
@ -838,8 +847,9 @@ CONTAINS
IF (.NOT. found) THEN
CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.FALSE., found=found, &
begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
END IF
END IF
CALL cif_get_real(parser, cell_angles(3))
cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
@ -1192,8 +1202,9 @@ CONTAINS
CALL parser_search_string(parser, "CRYST1", ignore_case=.FALSE., found=found, &
begin_line=.TRUE., search_from_begin_of_file=.TRUE.)
IF (.NOT. found) &
IF (.NOT. found) THEN
CPABORT("The line <CRYST1> was not found in PDB file! ")
END IF
periodic = 1
READ (parser%input_line, *, IOSTAT=ios) cryst, cell_lengths(:), cell_angles(:)

93
src/colvar_methods.F
View file

@ -618,11 +618,12 @@ CONTAINS
ALLOCATE (colvar%combine_cvs_param%variables(SIZE(my_par)))
colvar%combine_cvs_param%variables = my_par
! Check that the number of COLVAR provided is equal to the number of variables..
IF (SIZE(my_par) /= ncol) &
IF (SIZE(my_par) /= ncol) THEN
CALL cp_abort(__LOCATION__, &
"Number of defined COLVAR for COMBINE_COLVAR is different from the "// &
"number of variables! It is not possible to define COLVARs in a COMBINE_COLVAR "// &
"and avoid their usage in the combininig function!")
END IF
! Parameters
ALLOCATE (colvar%combine_cvs_param%c_parameters(0))
CALL section_vals_val_get(combine_section, "PARAMETERS", n_rep_val=ncol)
@ -726,8 +727,9 @@ CONTAINS
! Read the specification of the two planes
plane_sections => section_vals_get_subs_vals(plane_plane_angle_section, "PLANE")
CALL section_vals_get(plane_sections, n_repetition=n_var)
IF (n_var /= 2) &
IF (n_var /= 2) THEN
CPABORT("PLANE_PLANE_ANGLE Colvar section: Two PLANE sections must be provided!")
END IF
! Plane 1
CALL section_vals_val_get(plane_sections, "DEF_TYPE", i_rep_section=1, &
i_val=colvar%plane_plane_angle_param%plane1%type_of_def)
@ -736,8 +738,9 @@ CONTAINS
r_vals=s1)
colvar%plane_plane_angle_param%plane1%normal_vec = s1
IF (PRESENT(cell)) THEN
IF (ASSOCIATED(cell)) &
IF (ASSOCIATED(cell)) THEN
CALL cell_transform_input_cartesian(cell, colvar%plane_plane_angle_param%plane1%normal_vec)
END IF
END IF
ELSE
CALL section_vals_val_get(plane_sections, "ATOMS", i_rep_section=1, &
@ -753,8 +756,9 @@ CONTAINS
r_vals=s1)
colvar%plane_plane_angle_param%plane2%normal_vec = s1
IF (PRESENT(cell)) THEN
IF (ASSOCIATED(cell)) &
IF (ASSOCIATED(cell)) THEN
CALL cell_transform_input_cartesian(cell, colvar%plane_plane_angle_param%plane2%normal_vec)
END IF
END IF
ELSE
CALL section_vals_val_get(plane_sections, "ATOMS", i_rep_section=2, &
@ -863,9 +867,10 @@ CONTAINS
weights(ndim + 1:ndim + SIZE(wei)) = wei
ndim = ndim + SIZE(wei)
END DO
IF (ndim /= colvar%rmsd_param%n_atoms) &
IF (ndim /= colvar%rmsd_param%n_atoms) THEN
CALL cp_abort(__LOCATION__, "CV RMSD: list of atoms and list of "// &
"weights need to contain same number of entries. ")
END IF
DO i = 1, ndim
ii = colvar%rmsd_param%i_rmsd(i)
colvar%rmsd_param%weights(ii) = weights(i)
@ -947,11 +952,13 @@ CONTAINS
i_val=colvar%ring_puckering_param%iq)
! test the validity of the parameters
ndim = colvar%ring_puckering_param%nring
IF (ndim <= 3) &
IF (ndim <= 3) THEN
CPABORT("CV Ring Puckering: Ring size has to be 4 or larger. ")
END IF
ii = colvar%ring_puckering_param%iq
IF (ABS(ii) == 1 .OR. ii < -(ndim - 1)/2 .OR. ii > ndim/2) &
IF (ABS(ii) == 1 .OR. ii < -(ndim - 1)/2 .OR. ii > ndim/2) THEN
CPABORT("CV Ring Puckering: Invalid coordinate number.")
END IF
ELSE IF (my_subsection(23)) THEN
! Minimum Distance
wrk_section => mindist_section
@ -1326,7 +1333,7 @@ CONTAINS
IF (colvar%ring_puckering_param%iq == 0) THEN
WRITE (iw, '( A,T40,A)') ' COLVARS| Ring Puckering >>> coordinate', &
' Total Puckering Amplitude'
ELSEIF (colvar%ring_puckering_param%iq > 0) THEN
ELSE IF (colvar%ring_puckering_param%iq > 0) THEN
WRITE (iw, '( A,T35,A,T57,I8)') ' COLVARS| Ring Puckering >>> coordinate', &
' Puckering Amplitude', &
colvar%ring_puckering_param%iq
@ -2189,11 +2196,12 @@ CONTAINS
CALL put_derivative(colvar, iatom, fi)
END DO
ELSE
IF (force_env%in_use /= use_mixed_force) &
IF (force_env%in_use /= use_mixed_force) THEN
CALL cp_abort(__LOCATION__, &
'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
' A combination of mixed force_eval energies has been requested as '// &
' collective variable, but the MIXED env is not in use! Aborting.')
END IF
CALL force_env_get(force_env, force_env_section=force_env_section)
mapping_section => section_vals_get_subs_vals(force_env_section, "MIXED%MAPPING")
NULLIFY (values, parameters, subsystems, particles, global_forces, map_index, glob_natoms)
@ -2683,8 +2691,8 @@ CONTAINS
ss = ss - NINT(ss)
xkj = MATMUL(cell%hmat, ss)
! evaluation of the angle..
a = SQRT(DOT_PRODUCT(xij, xij))
b = SQRT(DOT_PRODUCT(xkj, xkj))
a = NORM2(xij)
b = NORM2(xkj)
t0 = 1.0_dp/(a*b)
t1 = 1.0_dp/(a**3.0_dp*b)
t2 = 1.0_dp/(a*b**3.0_dp)
@ -2834,8 +2842,8 @@ CONTAINS
ss = ss - NINT(ss)
xkj = MATMUL(cell%hmat, ss)
! Evaluation of the angle..
a = SQRT(DOT_PRODUCT(xij, xij))
b = SQRT(DOT_PRODUCT(xkj, xkj))
a = NORM2(xij)
b = NORM2(xkj)
t0 = 1.0_dp/(a*b)
t1 = 1.0_dp/(a**3.0_dp*b)
t2 = 1.0_dp/(a*b**3.0_dp)
@ -3137,10 +3145,6 @@ CONTAINS
TYPE(particle_list_type), POINTER :: particles_i
TYPE(particle_type), DIMENSION(:), POINTER :: my_particles
! settings for numerical derivatives
!REAL(KIND=dp) :: ri_step, dx_bond_j, dy_bond_j, dz_bond_j
!INTEGER :: idel
n_atoms_to = colvar%qparm_param%n_atoms_to
n_atoms_from = colvar%qparm_param%n_atoms_from
rcut = colvar%qparm_param%rcut
@ -3159,10 +3163,6 @@ CONTAINS
CPASSERT(r1cut < rcut)
denominator_tolerance = 1.0E-8_dp
!ri_step=0.1
!DO idel=-50, 50
!ftmp(:) = 0.0_dp
qparm = 0.0_dp
inv_n_atoms_from = 1.0_dp/REAL(n_atoms_from, KIND=dp)
DO ii = 1, n_atoms_from
@ -3200,11 +3200,10 @@ CONTAINS
shift(:) = 0.0_dp
shift(idim) = 1.0_dp
xij_shift = MATMUL(cell%hmat, shift)
rij_shift = SQRT(DOT_PRODUCT(xij_shift, xij_shift))
rij_shift = NORM2(xij_shift)
ncells(idim) = FLOOR(rcut/rij_shift - 0.5)
END DO !idim
!IF (mm.eq.0) WRITE(*,'(A8,3I3,A3,I10)') "Ncells:", ncells, "J:", j
shift(1:3) = 0.0_dp
DO aa = -ncells(1), ncells(1)
DO bb = -ncells(2), ncells(2)
@ -3215,12 +3214,7 @@ CONTAINS
shift(2) = REAL(bb, KIND=dp)
shift(3) = REAL(cc, KIND=dp)
xij = MATMUL(cell%hmat, ss0(:) + shift(:))
rij = SQRT(DOT_PRODUCT(xij, xij))
!IF (rij > rcut) THEN
! IF (mm==0) WRITE(*,'(A8,4F10.5)') " --", shift, rij
!ELSE
! IF (mm==0) WRITE(*,'(A8,4F10.5)') " ++", shift, rij
!ENDIF
rij = NORM2(xij)
IF (rij > rcut) CYCLE
! update qlm
@ -3236,7 +3230,7 @@ CONTAINS
IF (i == j) CYCLE jloop
xij(:) = xpj(:) - xpi(:)
rij = SQRT(DOT_PRODUCT(xij, xij))
rij = NORM2(xij)
IF (rij > rcut) CYCLE jloop
! update qlm
@ -3251,11 +3245,6 @@ CONTAINS
! this factor is necessary if one whishes to sum over m=0,L
! instead of m=-L,+L. This is off now because it is cheap and safe
fact = 1.0_dp
!IF (ABS(mm) > 0) THEN
! fact = 2.0_dp
!ELSE
! fact = 1.0_dp
!ENDIF
IF (nbond < denominator_tolerance) THEN
CPWARN("QPARM: number of neighbors is very close to zero!")
@ -3274,7 +3263,6 @@ CONTAINS
END DO ! loop over m
pre_fac = (4.0_dp*pi)/(2.0_dp*l + 1)
!WRITE(*,'(A8,2F10.5)') " si = ", SQRT(pre_fac*ql)
qparm = qparm + SQRT(pre_fac*ql)
ftmp(:) = 0.5_dp*SQRT(pre_fac/ql)*d_ql_dxi(:)
! multiply by -1 because aparently we have to save the force, not the gradient
@ -3287,10 +3275,6 @@ CONTAINS
colvar%ss = qparm*inv_n_atoms_from
colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
!WRITE(*,'(A15,3E20.10)') "COLVAR+DER = ", ri_step*idel, colvar%ss, -ftmp(1)
!ENDDO ! numercal derivative
END SUBROUTINE qparm_colvar
! **************************************************************************************************
@ -3323,7 +3307,6 @@ CONTAINS
exp_fac, fi, plm, pre_fac, sqrt_c1
REAL(KIND=dp), DIMENSION(3) :: dcosTheta, dfi
!bond = 1.0_dp/(1.0_dp+EXP(alpha*(rij-rcut)))
! RZK: infinitely differentiable smooth cutoff function
! that is precisely 1.0 below r1cut and precisely 0.0 above rcut
IF (rij > rcut) THEN
@ -3366,18 +3349,10 @@ CONTAINS
sqrt_c1 = SQRT(((2*ll + 1)*fac(ll - ABS(mm)))/(4*pi*fac(ll + ABS(mm))))
pre_fac = bond*sqrt_c1
dylm = pre_fac*dplm
!WHY? IF (plm < 0.0_dp) THEN
!WHY? dylm = -pre_fac*dplm
!WHY? ELSE
!WHY? dylm = pre_fac*dplm
!WHY? ENDIF
re_qlm = re_qlm + pre_fac*plm*COS(mm*fi)
im_qlm = im_qlm + pre_fac*plm*SIN(mm*fi)
!WRITE(*,'(A8,2I4,F10.5)') " Qlm = ", mm, j, bond
!WRITE(*,'(A8,2I4,2F10.5)') " Qlm = ", mm, j, re_qlm, im_qlm
dcosTheta(:) = xij(:)*xij(3)/(rij**3)
dcosTheta(3) = dcosTheta(3) - 1.0_dp/rij
! use tangent half-angle formula to compute d_fi/d_xi
@ -4595,7 +4570,7 @@ CONTAINS
IF (colvar%reaction_path_param%dist_rmsd) THEN
CALL rpath_dist_rmsd(colvar, my_particles)
ELSEIF (colvar%reaction_path_param%rmsd) THEN
ELSE IF (colvar%reaction_path_param%rmsd) THEN
CALL rpath_rmsd(colvar, my_particles)
ELSE
CALL rpath_colvar(colvar, cell, my_particles)
@ -4948,7 +4923,7 @@ CONTAINS
IF (colvar%reaction_path_param%dist_rmsd) THEN
CALL dpath_dist_rmsd(colvar, my_particles)
ELSEIF (colvar%reaction_path_param%rmsd) THEN
ELSE IF (colvar%reaction_path_param%rmsd) THEN
CALL dpath_rmsd(colvar, my_particles)
ELSE
CALL dpath_colvar(colvar, cell, my_particles)
@ -5864,11 +5839,12 @@ CONTAINS
DO j = 1, natom
! Atom coordinates
CALL parser_get_next_line(parser, 1, at_end=my_end)
IF (my_end) &
IF (my_end) THEN
CALL cp_abort(__LOCATION__, &
"Number of lines in XYZ format not equal to the number of atoms."// &
" Error in XYZ format for COORD_A (CV rmsd). Very probably the"// &
" line with title is missing or is empty. Please check the XYZ file and rerun your job!")
END IF
READ (parser%input_line, *) dummy_char, rptr(1:3)
r_ref((j - 1)*3 + 1, i) = cp_unit_to_cp2k(rptr(1), "angstrom")
r_ref((j - 1)*3 + 2, i) = cp_unit_to_cp2k(rptr(2), "angstrom")
@ -5966,12 +5942,7 @@ CONTAINS
iamin = wcai(i)
END IF
END DO
! zero=0.0_dp
! CALL put_derivative(colvar, 1, zero)
! CALL put_derivative(colvar, 2,zero)
! CALL put_derivative(colvar, 3, zero)
! write(*,'(2(i0,1x),4(f16.8,1x))')idmin,iamin,wc(1)%WannierHamDiag(idmin),wc(1)%WannierHamDiag(iamin),dmin,amin
colvar%ss = wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
DEALLOCATE (wcai)
DEALLOCATE (wcdi)
@ -5988,7 +5959,7 @@ CONTAINS
s = MATMUL(cell%h_inv, rij)
s = s - NINT(s)
xv = MATMUL(cell%hmat, s)
distance = SQRT(DOT_PRODUCT(xv, xv))
distance = NORM2(xv)
END FUNCTION distance
END SUBROUTINE Wc_colvar
@ -6007,8 +5978,8 @@ CONTAINS
TYPE(cell_type), POINTER :: cell
TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys
TYPE(particle_type), DIMENSION(:), &
OPTIONAL, POINTER :: particles
TYPE(qs_environment_type), POINTER, OPTIONAL :: qs_env ! optional just because I am lazy... but I should get rid of it...
OPTIONAL, POINTER :: particles
TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
INTEGER :: Od, H, Oa
REAL(dp) :: rOd(3), rOa(3), rH(3), &
@ -6105,7 +6076,7 @@ CONTAINS
s = MATMUL(cell%h_inv, rij)
s = s - NINT(s)
xv = MATMUL(cell%hmat, s)
distance = SQRT(DOT_PRODUCT(xv, xv))
distance = NORM2(xv)
END FUNCTION distance
END SUBROUTINE HBP_colvar

6
src/common/array_sort.fypp
View file

@ -97,9 +97,9 @@
! Merge will be performed directly in arr. Need backup of first sublist.
tmp_arr(1:m) = arr(1:m)
tmp_idx(1:m) = indices(1:m)
i = 1; ! number of elemens consumed from 1st sublist
j = 1; ! number of elemens consumed from 2nd sublist
k = 1; ! number of elemens already merged
i = 1 ! number of elements consumed from 1st sublist
j = 1 ! number of elements consumed from 2nd sublist
k = 1 ! number of elements already merged
DO WHILE (i <= m .and. j <= size(arr) - m)
IF (${prefix}$_less_than(arr(m + j), tmp_arr(i))) THEN

9
src/common/cp_array_utils.F
View file

@ -152,8 +152,9 @@ CONTAINS
WRITE (unit=unit_nr, fmt="(',')", advance="no")
END IF
END DO
IF (SIZE(array) > 0) &
IF (SIZE(array) > 0) THEN
WRITE (unit=unit_nr, fmt=el_format, advance="no") array(SIZE(array))
END IF
ELSE
DO i = 1, SIZE(array) - 1
WRITE (unit=unit_nr, fmt=defaultFormat, advance="no") array(i)
@ -163,8 +164,9 @@ CONTAINS
WRITE (unit=unit_nr, fmt="(',')", advance="no")
END IF
END DO
IF (SIZE(array) > 0) &
IF (SIZE(array) > 0) THEN
WRITE (unit=unit_nr, fmt=defaultFormat, advance="no") array(SIZE(array))
END IF
END IF
WRITE (unit=unit_nr, fmt="(' )')")
call m_flush(unit_nr)
@ -289,8 +291,7 @@ CONTAINS
!> \note
!> the array should be ordered in growing order
! **************************************************************************************************
FUNCTION cp_1d_${nametype1}$_bsearch(array, el, l_index, u_index) &
result(res)
FUNCTION cp_1d_${nametype1}$_bsearch(array, el, l_index, u_index) result(res)
${type1}$, intent(in) :: array(:)
${type1}$, intent(in) :: el
INTEGER, INTENT(in), OPTIONAL :: l_index, u_index

10
src/common/cp_error_handling.F
View file

@ -63,8 +63,9 @@ CONTAINS
CALL delay_non_master() ! cleaner output if all ranks abort simultaneously
unit_nr = cp_logger_get_default_io_unit()
IF (unit_nr <= 0) &
unit_nr = default_output_unit ! fall back to stdout
IF (unit_nr <= 0) THEN
unit_nr = default_output_unit
END IF ! fall back to stdout
CALL print_abort_message(message, location, unit_nr)
CALL print_stack(unit_nr)
@ -125,8 +126,9 @@ CONTAINS
! we (ab)use the logger to determine the first MPI rank
unit_nr = cp_logger_get_default_io_unit()
IF (unit_nr <= 0) &
wait_time = wait_time + 1.0_dp ! rank-0 gets a head start of one second.
IF (unit_nr <= 0) THEN
wait_time = wait_time + 1.0_dp
END IF ! rank-0 gets a head start of one second.
!$ IF (omp_get_thread_num() /= 0) &
!$ wait_time = wait_time + 1.0_dp ! master threads gets another second

33
src/common/cp_log_handling.F
View file

@ -303,8 +303,9 @@ CONTAINS
logger%ref_count = 1
IF (PRESENT(template_logger)) THEN
IF (template_logger%ref_count < 1) &
IF (template_logger%ref_count < 1) THEN
CPABORT(routineP//" template_logger%ref_count<1")
END IF
logger%print_level = template_logger%print_level
logger%default_global_unit_nr = template_logger%default_global_unit_nr
logger%close_local_unit_on_dealloc = template_logger%close_local_unit_on_dealloc
@ -339,16 +340,19 @@ CONTAINS
logger%suffix = ""
END IF
IF (PRESENT(para_env)) logger%para_env => para_env
IF (.NOT. ASSOCIATED(logger%para_env)) &
IF (.NOT. ASSOCIATED(logger%para_env)) THEN
CPABORT(routineP//" para env not associated")
IF (.NOT. logger%para_env%is_valid()) &
END IF
IF (.NOT. logger%para_env%is_valid()) THEN
CPABORT(routineP//" para_env%ref_count<1")
END IF
CALL logger%para_env%retain()
IF (PRESENT(print_level)) logger%print_level = print_level
IF (PRESENT(default_global_unit_nr)) &
IF (PRESENT(default_global_unit_nr)) THEN
logger%default_global_unit_nr = default_global_unit_nr
END IF
IF (PRESENT(global_filename)) THEN
logger%global_filename = global_filename
logger%close_global_unit_on_dealloc = .TRUE.
@ -362,8 +366,9 @@ CONTAINS
END IF
END IF
IF (PRESENT(default_local_unit_nr)) &
IF (PRESENT(default_local_unit_nr)) THEN
logger%default_local_unit_nr = default_local_unit_nr
END IF
IF (PRESENT(local_filename)) THEN
logger%local_filename = local_filename
logger%close_local_unit_on_dealloc = .TRUE.
@ -407,8 +412,9 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'cp_logger_retain', &
routineP = moduleN//':'//routineN
IF (logger%ref_count < 1) &
IF (logger%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
logger%ref_count = logger%ref_count + 1
END SUBROUTINE cp_logger_retain
@ -426,8 +432,9 @@ CONTAINS
routineP = moduleN//':'//routineN
IF (ASSOCIATED(logger)) THEN
IF (logger%ref_count < 1) &
IF (logger%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
logger%ref_count = logger%ref_count - 1
IF (logger%ref_count == 0) THEN
IF (logger%close_global_unit_on_dealloc .AND. &
@ -476,8 +483,9 @@ CONTAINS
lggr => logger
IF (.NOT. ASSOCIATED(lggr)) lggr => cp_get_default_logger()
IF (lggr%ref_count < 1) &
IF (lggr%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
res = level >= lggr%print_level
END FUNCTION cp_logger_would_log
@ -547,8 +555,9 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'cp_logger_set_log_level', &
routineP = moduleN//':'//routineN
IF (logger%ref_count < 1) &
IF (logger%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
logger%print_level = level
END SUBROUTINE cp_logger_set_log_level
@ -584,8 +593,9 @@ CONTAINS
NULLIFY (lggr)
END IF
IF (.NOT. ASSOCIATED(lggr)) lggr => cp_get_default_logger()
IF (lggr%ref_count < 1) &
IF (lggr%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
IF (PRESENT(local)) loc = local
IF (PRESENT(skip_not_ionode)) skip = skip_not_ionode
@ -692,8 +702,9 @@ CONTAINS
lggr => logger
IF (.NOT. ASSOCIATED(lggr)) lggr => cp_get_default_logger()
IF (lggr%ref_count < 1) &
IF (lggr%ref_count < 1) THEN
CPABORT(routineP//" logger%ref_count<1")
END IF
IF (PRESENT(local)) loc = local
IF (loc) THEN
res = TRIM(root)//TRIM(lggr%suffix)//'_p'// &

12
src/common/cp_result_methods.F
View file

@ -178,14 +178,16 @@ CONTAINS
n_rep = nrep
END IF
IF (nrep <= 0) &
IF (nrep <= 0) THEN
CALL cp_abort(__LOCATION__, &
" Trying to access result ("//TRIM(description)//") which was never stored!")
END IF
DO i = 1, nlist
IF (TRIM(results%result_label(i)) == TRIM(description)) THEN
IF (results%result_value(i)%value%type_in_use /= result_type_real) &
IF (results%result_value(i)%value%type_in_use /= result_type_real) THEN
CPABORT("Attempt to retrieve a RESULT which is not a REAL!")
END IF
size_res = SIZE(results%result_value(i)%value%real_type)
EXIT
@ -251,14 +253,16 @@ CONTAINS
n_rep = nrep
END IF
IF (nrep <= 0) &
IF (nrep <= 0) THEN
CALL cp_abort(__LOCATION__, &
" Trying to access result ("//TRIM(description)//") which was never stored!")
END IF
DO i = 1, nlist
IF (TRIM(results%result_label(i)) == TRIM(description)) THEN
IF (results%result_value(i)%value%type_in_use /= result_type_real) &
IF (results%result_value(i)%value%type_in_use /= result_type_real) THEN
CPABORT("Attempt to retrieve a RESULT which is not a REAL!")
END IF
size_res = SIZE(results%result_value(i)%value%real_type)
EXIT

36
src/common/cp_units.F
View file

@ -650,10 +650,11 @@ CONTAINS
CPABORT("unknown electric field unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
CASE (cp_ukind_none)
IF (basic_unit /= cp_units_none) &
IF (basic_unit /= cp_units_none) THEN
CALL cp_abort(__LOCATION__, &
"if the kind of the unit is none also unit must be undefined,not:" &
//TRIM(cp_to_string(basic_unit)))
END IF
CASE default
CPABORT("unknown kind of unit:"//TRIM(cp_to_string(basic_kind)))
END SELECT
@ -679,10 +680,11 @@ CONTAINS
my_power = 1
IF (PRESENT(power)) my_power = power
IF (basic_unit == cp_units_none .AND. basic_kind /= cp_ukind_undef) THEN
IF (basic_kind /= cp_units_none) &
IF (basic_kind /= cp_units_none) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(cp_to_string(basic_unit)))
END IF
END IF
SELECT CASE (basic_kind)
CASE (cp_ukind_undef)
@ -862,9 +864,10 @@ CONTAINS
IF (accept_undefined) my_accept_undefined = accept_undefined
IF (PRESENT(power)) my_power = power
IF (basic_unit == cp_units_none) THEN
IF (.NOT. my_accept_undefined .AND. basic_kind == cp_units_none) &
IF (.NOT. my_accept_undefined .AND. basic_kind == cp_units_none) THEN
CALL cp_abort(__LOCATION__, "unit not yet fully specified, unit of kind "// &
TRIM(cp_to_string(basic_kind)))
END IF
END IF
SELECT CASE (basic_kind)
CASE (cp_ukind_undef)
@ -900,10 +903,11 @@ CONTAINS
res = "K_e"
CASE (cp_units_none)
res = "energy"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown energy unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -931,10 +935,11 @@ CONTAINS
res = "au_temp"
CASE (cp_units_none)
res = "temperature"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown temperature unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -956,10 +961,11 @@ CONTAINS
res = "au_p"
CASE (cp_units_none)
res = "pressure"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown pressure unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -971,10 +977,11 @@ CONTAINS
res = "deg"
CASE (cp_units_none)
res = "angle"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown angle unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -992,10 +999,11 @@ CONTAINS
res = "wavenumber_t"
CASE (cp_units_none)
res = "time"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown time unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -1009,10 +1017,11 @@ CONTAINS
res = "m_e"
CASE (cp_units_none)
res = "mass"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown mass unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -1024,10 +1033,11 @@ CONTAINS
res = "au_pot"
CASE (cp_units_none)
res = "potential"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown potential unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -1041,10 +1051,11 @@ CONTAINS
res = "au_f"
CASE (cp_units_none)
res = "force"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown potential unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT
@ -1060,10 +1071,11 @@ CONTAINS
res = "au_efield"
CASE (cp_units_none)
res = "electric field"
IF (.NOT. my_accept_undefined) &
IF (.NOT. my_accept_undefined) THEN
CALL cp_abort(__LOCATION__, &
"unit not yet fully specified, unit of kind "// &
TRIM(res))
END IF
CASE default
CPABORT("unknown efield unit:"//TRIM(cp_to_string(basic_unit)))
END SELECT

3
src/common/distribution_1d_types.F
View file

@ -104,8 +104,9 @@ CONTAINS
CALL para_env%retain()
distribution_1d%listbased_distribution = .FALSE.
IF (PRESENT(listbased_distribution)) &
IF (PRESENT(listbased_distribution)) THEN
distribution_1d%listbased_distribution = listbased_distribution
END IF
ALLOCATE (distribution_1d%n_el(my_n_lists), distribution_1d%list(my_n_lists))

34
src/common/fparser.F
View file

@ -407,10 +407,10 @@ CONTAINS
DO j = b, e
IF (Func(j:j) == '(') THEN
ParCnt = ParCnt + 1
ELSEIF (Func(j:j) == ')') THEN
ELSE IF (Func(j:j) == ')') THEN
ParCnt = ParCnt - 1
IF (ParCnt == 0) EXIT
ELSEIF (ParCnt == 1 .AND. Func(j:j) == ',') THEN
ELSE IF (ParCnt == 1 .AND. Func(j:j) == ',') THEN
ArgPos = j
ArgCnt = ArgCnt + 1
END IF
@ -748,7 +748,7 @@ CONTAINS
DO j = b + 1, e - 1
IF (F(j:j) == '(') THEN
k = k + 1
ELSEIF (F(j:j) == ')') THEN
ELSE IF (F(j:j) == ')') THEN
k = k - 1
END IF
IF (k < 0) EXIT
@ -792,11 +792,11 @@ CONTAINS
! WRITE(*,*)'1. F(b:e) = "+..."'
CALL CompileSubstr(i, F, b + 1, e, Var)
RETURN
ELSEIF (CompletelyEnclosed(F, b, e)) THEN ! Case 2: F(b:e) = '(...)'
ELSE IF (CompletelyEnclosed(F, b, e)) THEN ! Case 2: F(b:e) = '(...)'
! WRITE(*,*)'2. F(b:e) = "(...)"'
CALL CompileSubstr(i, F, b + 1, e - 1, Var)
RETURN
ELSEIF (SCAN(F(b:b), calpha) > 0) THEN
ELSE IF (SCAN(F(b:b), calpha) > 0) THEN
n = MathFunctionIndex(F(b:e))
IF (n > 0) THEN
b2 = b + INDEX(F(b:e), '(') - 1
@ -807,13 +807,13 @@ CONTAINS
RETURN
END IF
END IF
ELSEIF (F(b:b) == '-') THEN
ELSE IF (F(b:b) == '-') THEN
IF (CompletelyEnclosed(F, b + 1, e)) THEN ! Case 4: F(b:e) = '-(...)'
! WRITE(*,*)'4. F(b:e) = "-(...)"'
CALL CompileSubstr(i, F, b + 2, e - 1, Var)
CALL AddCompiledByte(i, cNeg)
RETURN
ELSEIF (SCAN(F(b + 1:b + 1), calpha) > 0) THEN
ELSE IF (SCAN(F(b + 1:b + 1), calpha) > 0) THEN
n = MathFunctionIndex(F(b + 1:e))
IF (n > 0) THEN
b2 = b + INDEX(F(b + 1:e), '(')
@ -835,7 +835,7 @@ CONTAINS
DO j = e, b, -1
IF (F(j:j) == ')') THEN
k = k + 1
ELSEIF (F(j:j) == '(') THEN
ELSE IF (F(j:j) == '(') THEN
k = k - 1
END IF
IF (k == 0 .AND. F(j:j) == Ops(io) .AND. IsBinaryOp(j, F)) THEN
@ -900,9 +900,9 @@ CONTAINS
DO j = b, e
IF (F(j:j) == '(') THEN
ParCnt = ParCnt + 1
ELSEIF (F(j:j) == ')') THEN
ELSE IF (F(j:j) == ')') THEN
ParCnt = ParCnt - 1
ELSEIF (ParCnt == 0 .AND. F(j:j) == ',') THEN
ELSE IF (ParCnt == 0 .AND. F(j:j) == ',') THEN
CALL CompileSubstr(i, F, b2, j - 1, Var)
b2 = j + 1
END IF
@ -938,17 +938,17 @@ CONTAINS
IF (F(j:j) == '+' .OR. F(j:j) == '-') THEN ! Plus or minus sign:
IF (j == 1) THEN ! - leading unary operator ?
res = .FALSE.
ELSEIF (SCAN(F(j - 1:j - 1), '+-*/^(,') > 0) THEN ! - other unary operator ?
ELSE IF (SCAN(F(j - 1:j - 1), '+-*/^(,') > 0) THEN ! - other unary operator ?
res = .FALSE.
ELSEIF (SCAN(F(j + 1:j + 1), '0123456789') > 0 .AND. & ! - in exponent of real number ?
SCAN(F(j - 1:j - 1), 'eEdD') > 0) THEN
ELSE IF (SCAN(F(j + 1:j + 1), '0123456789') > 0 .AND. & ! - in exponent of real number ?
SCAN(F(j - 1:j - 1), 'eEdD') > 0) THEN
Dflag = .FALSE.; Pflag = .FALSE.
k = j - 1
DO WHILE (k > 1) ! step to the left in mantissa
k = k - 1
IF (SCAN(F(k:k), '0123456789') > 0) THEN
Dflag = .TRUE.
ELSEIF (F(k:k) == '.') THEN
ELSE IF (F(k:k) == '.') THEN
IF (Pflag) THEN
EXIT ! * EXIT: 2nd appearance of '.'
ELSE
@ -1011,7 +1011,7 @@ CONTAINS
CASE ('+', '-') ! Permitted only
IF (Bflag) THEN
InMan = .TRUE.; Bflag = .FALSE. ! - at beginning of mantissa
ELSEIF (Eflag) THEN
ELSE IF (Eflag) THEN
InExp = .TRUE.; Eflag = .FALSE. ! - at beginning of exponent
ELSE
EXIT ! - otherwise STOP
@ -1019,7 +1019,7 @@ CONTAINS
CASE ('0':'9') ! Mark
IF (Bflag) THEN
InMan = .TRUE.; Bflag = .FALSE. ! - beginning of mantissa
ELSEIF (Eflag) THEN
ELSE IF (Eflag) THEN
InExp = .TRUE.; Eflag = .FALSE. ! - beginning of exponent
END IF
IF (InMan) DInMan = .TRUE. ! Mantissa contains digit
@ -1028,7 +1028,7 @@ CONTAINS
IF (Bflag) THEN
Pflag = .TRUE. ! - mark 1st appearance of '.'
InMan = .TRUE.; Bflag = .FALSE. ! mark beginning of mantissa
ELSEIF (InMan .AND. .NOT. Pflag) THEN
ELSE IF (InMan .AND. .NOT. Pflag) THEN
Pflag = .TRUE. ! - mark 1st appearance of '.'
ELSE
EXIT ! - otherwise STOP

2
src/common/gfun.F
View file

@ -62,7 +62,7 @@ CONTAINS
IF (t <= 12.0_dp) THEN
! downward recursion
g(nmax) = gfun_taylor(nmax, t);
g(nmax) = gfun_taylor(nmax, t)
DO i = nmax, 1, -1
g(i - 1) = (1.0_dp - 2.0_dp*t*g(i))/(2.0_dp*i - 1.0_dp)
END DO

6
src/common/hash_map.fypp
View file

@ -81,11 +81,13 @@
initial_capacity_ = 11
END IF
IF (initial_capacity_ < 1) &
IF (initial_capacity_ < 1) THEN
CPABORT("initial_capacity < 1")
END IF
IF (ASSOCIATED(hash_map%buckets)) &
IF (ASSOCIATED(hash_map%buckets)) THEN
CPABORT("hash map is already initialized.")
END IF
ALLOCATE (hash_map%buckets(initial_capacity_))
hash_map%size = 0

114
src/common/list.fypp
View file

@ -87,15 +87,18 @@
initial_capacity_ = 11
If (PRESENT(initial_capacity)) initial_capacity_ = initial_capacity
IF (initial_capacity_ < 0) &
IF (initial_capacity_ < 0) THEN
CPABORT("list_${valuetype}$_create: initial_capacity < 0")
END IF
IF (ASSOCIATED(list%arr)) &
IF (ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_create: list is already initialized.")
END IF
ALLOCATE (list%arr(initial_capacity_), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("list_${valuetype}$_init: allocation failed")
END IF
list%size = 0
END SUBROUTINE list_${valuetype}$_init
@ -112,8 +115,9 @@
SUBROUTINE list_${valuetype}$_destroy(list)
TYPE(list_${valuetype}$_type), intent(inout) :: list
INTEGER :: i
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_destroy: list is not initialized.")
END IF
do i = 1, list%size
deallocate (list%arr(i)%p)
@ -137,12 +141,15 @@
TYPE(list_${valuetype}$_type), intent(inout) :: list
${valuetype_in}$, intent(in) :: value
INTEGER, intent(in) :: pos
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_set: list is not initialized.")
IF (pos < 1) &
END IF
IF (pos < 1) THEN
CPABORT("list_${valuetype}$_set: pos < 1")
IF (pos > list%size) &
END IF
IF (pos > list%size) THEN
CPABORT("list_${valuetype}$_set: pos > size")
END IF
list%arr(pos)%p%value ${value_assign}$value
END SUBROUTINE list_${valuetype}$_set
@ -159,15 +166,18 @@
${valuetype_in}$, intent(in) :: value
INTEGER :: stat
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_push: list is not initialized.")
if (list%size == size(list%arr)) &
call change_capacity_${valuetype}$ (list, 2*size(list%arr) + 1)
END IF
IF (list%size == size(list%arr)) THEN
CALL change_capacity_${valuetype}$ (list, 2*size(list%arr) + 1)
END IF
list%size = list%size + 1
ALLOCATE (list%arr(list%size)%p, stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("list_${valuetype}$_push: allocation failed")
END IF
list%arr(list%size)%p%value ${value_assign}$value
END SUBROUTINE list_${valuetype}$_push
@ -187,15 +197,19 @@
INTEGER, intent(in) :: pos
INTEGER :: i, stat
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_insert: list is not initialized.")
IF (pos < 1) &
END IF
IF (pos < 1) THEN
CPABORT("list_${valuetype}$_insert: pos < 1")
IF (pos > list%size + 1) &
END IF
IF (pos > list%size + 1) THEN
CPABORT("list_${valuetype}$_insert: pos > size+1")
END IF
if (list%size == size(list%arr)) &
if (list%size == size(list%arr)) THEN
call change_capacity_${valuetype}$ (list, 2*size(list%arr) + 1)
END IF
list%size = list%size + 1
do i = list%size, pos + 1, -1
@ -203,8 +217,9 @@
end do
ALLOCATE (list%arr(pos)%p, stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("list_${valuetype}$_insert: allocation failed.")
END IF
list%arr(pos)%p%value ${value_assign}$value
END SUBROUTINE list_${valuetype}$_insert
@ -221,10 +236,12 @@
TYPE(list_${valuetype}$_type), intent(inout) :: list
${valuetype_out}$ :: value
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_peek: list is not initialized.")
IF (list%size < 1) &
END IF
IF (list%size < 1) THEN
CPABORT("list_${valuetype}$_peek: list is empty.")
END IF
value ${value_assign}$list%arr(list%size)%p%value
END FUNCTION list_${valuetype}$_peek
@ -245,10 +262,12 @@
TYPE(list_${valuetype}$_type), intent(inout) :: list
${valuetype_out}$ :: value
IF (.not. ASSOCIATED(list%arr)) &
IF (.NOT. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_pop: list is not initialized.")
IF (list%size < 1) &
END IF
IF (list%size < 1) THEN
CPABORT("list_${valuetype}$_pop: list is empty.")
END IF
value ${value_assign}$list%arr(list%size)%p%value
deallocate (list%arr(list%size)%p)
@ -266,12 +285,13 @@
TYPE(list_${valuetype}$_type), intent(inout) :: list
INTEGER :: i
IF (.not. ASSOCIATED(list%arr)) &
IF (.not. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_clear: list is not initialized.")
END IF
do i = 1, list%size
DO i = 1, list%size
deallocate (list%arr(i)%p)
end do
END DO
list%size = 0
END SUBROUTINE list_${valuetype}$_clear
@ -290,12 +310,15 @@
INTEGER, intent(in) :: pos
${valuetype_out}$ :: value
IF (.not. ASSOCIATED(list%arr)) &
IF (.NOT. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_get: list is not initialized.")
IF (pos < 1) &
END IF
IF (pos < 1) THEN
CPABORT("list_${valuetype}$_get: pos < 1")
IF (pos > list%size) &
END IF
IF (pos > list%size) THEN
CPABORT("list_${valuetype}$_get: pos > size")
END IF
value ${value_assign}$list%arr(pos)%p%value
@ -314,17 +337,20 @@
INTEGER, intent(in) :: pos
INTEGER :: i
IF (.not. ASSOCIATED(list%arr)) &
IF (.NOT. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_del: list is not initialized.")
IF (pos < 1) &
END IF
IF (pos < 1) THEN
CPABORT("list_${valuetype}$_det: pos < 1")
IF (pos > list%size) &
END IF
IF (pos > list%size) THEN
CPABORT("list_${valuetype}$_det: pos > size")
END IF
deallocate (list%arr(pos)%p)
do i = pos, list%size - 1
DO i = pos, list%size - 1
list%arr(i)%p => list%arr(i + 1)%p
end do
END DO
list%size = list%size - 1
@ -342,8 +368,9 @@
TYPE(list_${valuetype}$_type), intent(in) :: list
INTEGER :: size
IF (.not. ASSOCIATED(list%arr)) &
IF (.NOT. ASSOCIATED(list%arr)) THEN
CPABORT("list_${valuetype}$_size: list is not initialized.")
END IF
size = list%size
END FUNCTION list_${valuetype}$_size
@ -363,29 +390,34 @@
TYPE(private_item_p_type_${valuetype}$), DIMENSION(:), POINTER :: old_arr
new_cap = new_capacity
IF (new_cap < 0) &
IF (new_cap < 0) THEN
CPABORT("list_${valuetype}$_change_capacity: new_capacity < 0")
IF (new_cap < list%size) &
END IF
IF (new_cap < list%size) THEN
CPABORT("list_${valuetype}$_change_capacity: new_capacity < size")
END IF
IF (new_cap > HUGE(i)) THEN
IF (size(list%arr) == HUGE(i)) &
IF (size(list%arr) == HUGE(i)) THEN
CPABORT("list_${valuetype}$_change_capacity: list has reached integer limit.")
END IF
new_cap = HUGE(i) ! grow as far as possible
END IF
old_arr => list%arr
allocate (list%arr(new_cap), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("list_${valuetype}$_change_capacity: allocation failed")
END IF
do i = 1, list%size
DO i = 1, list%size
allocate (list%arr(i)%p, stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("list_${valuetype}$_change_capacity: allocation failed")
END IF
list%arr(i)%p%value ${value_assign}$old_arr(i)%p%value
deallocate (old_arr(i)%p)
end do
deallocate (old_arr)
DEALLOCATE (old_arr(i)%p)
END DO
DEALLOCATE (old_arr)
END SUBROUTINE change_capacity_${valuetype}$
#:enddef

37
src/common/mathlib.F
View file

@ -187,8 +187,8 @@ CONTAINS
REAL(KIND=dp) :: length_of_a, length_of_b
REAL(KIND=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
length_of_a = SQRT(DOT_PRODUCT(a, a))
length_of_b = SQRT(DOT_PRODUCT(b, b))
length_of_a = NORM2(a)
length_of_b = NORM2(b)
IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
a_norm(:) = a(:)/length_of_a
@ -997,8 +997,9 @@ CONTAINS
! set singular values that are too small to zero
DO i = 1, n
IF (sig(i) > rskip*MAXVAL(sig)) THEN
IF (PRESENT(determinant)) &
IF (PRESENT(determinant)) THEN
determinant = determinant*sig(i)
END IF
sig_plus(i, i) = 1._dp/sig(i)
ELSE
sig_plus(i, i) = 0.0_dp
@ -1975,12 +1976,14 @@ CONTAINS
IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
A_trans == 'T' .OR. A_trans == 't' .OR. &
A_trans == 'C' .OR. A_trans == 'c')) &
A_trans == 'C' .OR. A_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 1 (A).")
END IF
IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
B_trans == 'T' .OR. B_trans == 't' .OR. &
B_trans == 'C' .OR. B_trans == 'c')) &
B_trans == 'C' .OR. B_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 2 (B).")
END IF
CALL timeset(routineN, handle)
@ -2018,12 +2021,14 @@ CONTAINS
IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
A_trans == 'T' .OR. A_trans == 't' .OR. &
A_trans == 'C' .OR. A_trans == 'c')) &
A_trans == 'C' .OR. A_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 1 (A).")
END IF
IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
B_trans == 'T' .OR. B_trans == 't' .OR. &
B_trans == 'C' .OR. B_trans == 'c')) &
B_trans == 'C' .OR. B_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 2 (B).")
END IF
CALL timeset(routineN, handle)
@ -2068,16 +2073,19 @@ CONTAINS
IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
A_trans == 'T' .OR. A_trans == 't' .OR. &
A_trans == 'C' .OR. A_trans == 'c')) &
A_trans == 'C' .OR. A_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 1 (A).")
END IF
IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
B_trans == 'T' .OR. B_trans == 't' .OR. &
B_trans == 'C' .OR. B_trans == 'c')) &
B_trans == 'C' .OR. B_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 2 (B).")
END IF
IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
C_trans == 'T' .OR. C_trans == 't' .OR. &
C_trans == 'C' .OR. C_trans == 'c')) &
C_trans == 'C' .OR. C_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 3 (C).")
END IF
CALL timeset(routineN, handle)
@ -2127,16 +2135,19 @@ CONTAINS
IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
A_trans == 'T' .OR. A_trans == 't' .OR. &
A_trans == 'C' .OR. A_trans == 'c')) &
A_trans == 'C' .OR. A_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 1 (A).")
END IF
IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
B_trans == 'T' .OR. B_trans == 't' .OR. &
B_trans == 'C' .OR. B_trans == 'c')) &
B_trans == 'C' .OR. B_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 2 (B).")
END IF
IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
C_trans == 'T' .OR. C_trans == 't' .OR. &
C_trans == 'C' .OR. C_trans == 'c')) &
C_trans == 'C' .OR. C_trans == 'c')) THEN
CPABORT("Unknown transpose character for array 3 (C).")
END IF
CALL timeset(routineN, handle)

27
src/common/memory_utilities_unittest.F
View file

@ -32,11 +32,13 @@ CONTAINS
CALL reallocate(real_arr, 1, 20)
IF (.NOT. ALL(real_arr(1:10) == [(idx, idx=1, 10)])) &
IF (.NOT. ALL(real_arr(1:10) == [(idx, idx=1, 10)])) THEN
ERROR STOP "check_real_rank1_allocated: reallocating changed the initial values"
END IF
IF (.NOT. ALL(real_arr(11:20) == 0.)) &
IF (.NOT. ALL(real_arr(11:20) == 0.)) THEN
ERROR STOP "check_real_rank1_allocated: reallocation failed to initialise new values with 0."
END IF
DEALLOCATE (real_arr)
@ -53,8 +55,9 @@ CONTAINS
CALL reallocate(real_arr, 1, 20)
IF (.NOT. ALL(real_arr(1:20) == 0.)) &
IF (.NOT. ALL(real_arr(1:20) == 0.)) THEN
ERROR STOP "check_real_rank1_unallocated: reallocation failed to initialise new values with 0."
END IF
DEALLOCATE (real_arr)
@ -73,11 +76,13 @@ CONTAINS
CALL reallocate(real_arr, 1, 10, 1, 5)
IF (.NOT. (ALL(real_arr(1:5, 1) == [(idx, idx=1, 5)]) .AND. ALL(real_arr(1:5, 2) == [(idx, idx=6, 10)]))) &
IF (.NOT. (ALL(real_arr(1:5, 1) == [(idx, idx=1, 5)]) .AND. ALL(real_arr(1:5, 2) == [(idx, idx=6, 10)]))) THEN
ERROR STOP "check_real_rank2_allocated: reallocating changed the initial values"
END IF
IF (.NOT. (ALL(real_arr(6:10, 1:2) == 0.) .AND. ALL(real_arr(1:10, 3:5) == 0.))) &
IF (.NOT. (ALL(real_arr(6:10, 1:2) == 0.) .AND. ALL(real_arr(1:10, 3:5) == 0.))) THEN
ERROR STOP "check_real_rank2_allocated: reallocation failed to initialise new values with 0."
END IF
DEALLOCATE (real_arr)
@ -94,8 +99,9 @@ CONTAINS
CALL reallocate(real_arr, 1, 10, 1, 5)
IF (.NOT. ALL(real_arr(1:10, 1:5) == 0.)) &
IF (.NOT. ALL(real_arr(1:10, 1:5) == 0.)) THEN
ERROR STOP "check_real_rank2_unallocated: reallocation failed to initialise new values with 0."
END IF
DEALLOCATE (real_arr)
@ -114,11 +120,13 @@ CONTAINS
CALL reallocate(str_arr, 1, 20)
IF (.NOT. ALL(str_arr(1:10) == [("hello, there", idx=1, 10)])) &
IF (.NOT. ALL(str_arr(1:10) == [("hello, there", idx=1, 10)])) THEN
ERROR STOP "check_string_rank1_allocated: reallocating changed the initial values"
END IF
IF (.NOT. ALL(str_arr(11:20) == "")) &
IF (.NOT. ALL(str_arr(11:20) == "")) THEN
ERROR STOP "check_string_rank1_allocated: reallocation failed to initialise new values with ''."
END IF
DEALLOCATE (str_arr)
@ -135,8 +143,9 @@ CONTAINS
CALL reallocate(str_arr, 1, 20)
IF (.NOT. ALL(str_arr(1:20) == "")) &
IF (.NOT. ALL(str_arr(1:20) == "")) THEN
ERROR STOP "check_string_rank1_allocated: reallocation failed to initialise new values with ''."
END IF
DEALLOCATE (str_arr)

12
src/common/parallel_rng_types.F
View file

@ -512,8 +512,9 @@ CONTAINS
LOGICAL, INTENT(IN), OPTIONAL :: antithetic, extended_precision
TYPE(rng_stream_type) :: rng_stream
IF (LEN_TRIM(name) > rng_name_length) &
IF (LEN_TRIM(name) > rng_name_length) THEN
CPABORT("given random number generator name is too long")
END IF
rng_stream%name = TRIM(name)
@ -630,14 +631,16 @@ CONTAINS
LOGICAL, INTENT(OUT), OPTIONAL :: buffer_filled
IF (PRESENT(name)) name = self%name
IF (PRESENT(distribution_type)) &
IF (PRESENT(distribution_type)) THEN
distribution_type = self%distribution_type
END IF
IF (PRESENT(bg)) bg = self%bg
IF (PRESENT(cg)) cg = self%cg
IF (PRESENT(ig)) ig = self%ig
IF (PRESENT(antithetic)) antithetic = self%antithetic
IF (PRESENT(extended_precision)) &
IF (PRESENT(extended_precision)) THEN
extended_precision = self%extended_precision
END IF
IF (PRESENT(buffer)) buffer = self%buffer
IF (PRESENT(buffer_filled)) buffer_filled = self%buffer_filled
END SUBROUTINE get
@ -1131,8 +1134,9 @@ CONTAINS
my_write_all = .FALSE.
IF (PRESENT(write_all)) &
IF (PRESENT(write_all)) THEN
my_write_all = write_all
END IF
WRITE (UNIT=output_unit, FMT="(/,T2,A,/)") &
"Random number stream <"//TRIM(self%name)//">:"

30
src/common/parallel_rng_types_unittest.F
View file

@ -32,14 +32,16 @@ PROGRAM parallel_rng_types_TEST
nsamples = 1000
nargs = command_argument_count()
IF (nargs > 1) &
IF (nargs > 1) then
ERROR STOP "Usage: parallel_rng_types_TEST [<int:nsamples>]"
end if
IF (nargs == 1) THEN
CALL get_command_argument(1, arg)
READ (arg, *, iostat=stat) nsamples
IF (stat /= 0) &
IF (stat /= 0) then
ERROR STOP "Usage: parallel_rng_types_TEST [<int:nsamples>]"
end if
END IF
CALL mp_world_init(mpi_comm)
@ -60,8 +62,9 @@ PROGRAM parallel_rng_types_TEST
distribution_type=UNIFORM, &
extended_precision=.TRUE.)
IF (ionode) &
IF (ionode) then
CALL rng_stream%write(default_output_unit)
end if
tmax = -HUGE(0.0_dp)
tmin = +HUGE(0.0_dp)
@ -94,8 +97,9 @@ PROGRAM parallel_rng_types_TEST
distribution_type=GAUSSIAN, &
extended_precision=.TRUE.)
IF (ionode) &
IF (ionode) then
CALL rng_stream%write(default_output_unit)
end if
tmax = -HUGE(0.0_dp)
tmin = +HUGE(0.0_dp)
@ -162,8 +166,9 @@ CONTAINS
CALL rng_stream%get(ig=ig, cg=cg, bg=bg, name=name)
IF (ANY(ig /= ig_orig) .OR. ANY(cg /= cg_orig) .OR. ANY(bg /= bg_orig) &
.OR. (name /= name_orig)) &
.OR. (name /= name_orig)) then
ERROR STOP "Stream dump and load roundtrip failed"
end if
WRITE (UNIT=default_output_unit, FMT="(T4,A)") &
"Roundtrip successful"
@ -185,8 +190,9 @@ CONTAINS
WRITE (UNIT=default_output_unit, FMT="(T4,A10,A433)") &
"GENERATED:", rng_record
IF (rng_record /= serialized_string) &
IF (rng_record /= serialized_string) then
ERROR STOP "Serialized record does not match the expected output"
end if
WRITE (UNIT=default_output_unit, FMT="(T4,A)") &
"Serialized record matches the expected output"
@ -214,19 +220,22 @@ CONTAINS
arr = orig
CALL rng_stream%shuffle(arr)
IF (ALL(arr == orig)) &
IF (ALL(arr == orig)) then
ERROR STOP "shuffle failed: array was left untouched"
end if
WRITE (UNIT=default_output_unit, FMT="(A)", ADVANCE="no") "."
IF (ANY(arr /= orig(arr))) &
IF (ANY(arr /= orig(arr))) then
ERROR STOP "shuffle failed: the shuffled original is not the shuffled original"
end if
WRITE (UNIT=default_output_unit, FMT="(A)", ADVANCE="no") "."
! sort and compare to orig
mask = .TRUE.
DO idx = 1, size(orig)
IF (MINVAL(arr, mask) /= orig(idx)) &
IF (MINVAL(arr, mask) /= orig(idx)) then
ERROR STOP "shuffle failed: there is at least one unknown index"
end if
mask(MINLOC(arr, mask)) = .FALSE.
END DO
WRITE (UNIT=default_output_unit, FMT="(A)", ADVANCE="no") "."
@ -235,8 +244,9 @@ CONTAINS
CALL rng_stream%reset()
CALL rng_stream%shuffle(arr2)
IF (ANY(arr2 /= arr)) &
IF (ANY(arr2 /= arr)) then
ERROR STOP "shuffle failed: array was shuffled differently with same rng state"
end if
WRITE (UNIT=default_output_unit, FMT="(A)", ADVANCE="no") "."
WRITE (UNIT=default_output_unit, FMT="(T4,A)") &

12
src/common/reference_manager.F
View file

@ -332,14 +332,18 @@ CONTAINS
WRITE (unit, '(T3,A)') '<SOURCE>'//thebib(i)%ref%source//'</SOURCE>'
! DOI, volume, pages, year, month.
IF (ALLOCATED(thebib(i)%ref%doi)) &
IF (ALLOCATED(thebib(i)%ref%doi)) THEN
WRITE (unit, '(T3,A)') '<DOI>'//TRIM(substitute_special_xml_tokens(thebib(i)%ref%doi))//'</DOI>'
IF (ALLOCATED(thebib(i)%ref%volume)) &
END IF
IF (ALLOCATED(thebib(i)%ref%volume)) THEN
WRITE (unit, '(T3,A)') '<VOLUME>'//thebib(i)%ref%volume//'</VOLUME>'
IF (ALLOCATED(thebib(i)%ref%pages)) &
END IF
IF (ALLOCATED(thebib(i)%ref%pages)) THEN
WRITE (unit, '(T3,A)') '<PAGES>'//thebib(i)%ref%pages//'</PAGES>'
IF (thebib(i)%ref%year > 0) &
END IF
IF (thebib(i)%ref%year > 0) THEN
WRITE (unit, '(T3,A,I4.4,A)') '<YEAR>', thebib(i)%ref%year, '</YEAR>'
END IF
WRITE (unit, '(T2,A)') '</REFERENCE>'
END DO

40
src/common/spherical_harmonics.F
View file

@ -1805,30 +1805,42 @@ CONTAINS
REAL(KIND=dp) :: sumk, t
! Check validity of the input parameters
IF (j1 < 0.0_dp) &
IF (j1 < 0.0_dp) THEN
CPABORT("The angular momentum quantum number j1 has to be nonnegative")
IF (.NOT. (is_integer(j1) .OR. is_integer(2.0_dp*j1))) &
END IF
IF (.NOT. (is_integer(j1) .OR. is_integer(2.0_dp*j1))) THEN
CPABORT("The angular momentum quantum number j1 has to be integer or half-integer")
IF (j2 < 0.0_dp) &
END IF
IF (j2 < 0.0_dp) THEN
CPABORT("The angular momentum quantum number j2 has to be nonnegative")
IF (.NOT. (is_integer(j2) .OR. is_integer(2.0_dp*j2))) &
END IF
IF (.NOT. (is_integer(j2) .OR. is_integer(2.0_dp*j2))) THEN
CPABORT("The angular momentum quantum number j2 has to be integer or half-integer")
IF (J < 0.0_dp) &
END IF
IF (J < 0.0_dp) THEN
CPABORT("The angular momentum quantum number J has to be nonnegative")
IF (.NOT. (is_integer(J) .OR. is_integer(2.0_dp*J))) &
END IF
IF (.NOT. (is_integer(J) .OR. is_integer(2.0_dp*J))) THEN
CPABORT("The angular momentum quantum number J has to be integer or half-integer")
IF ((ABS(m1) - j1) > EPSILON(m1)) &
END IF
IF ((ABS(m1) - j1) > EPSILON(m1)) THEN
CPABORT("The angular momentum quantum number m1 has to satisfy -j1 <= m1 <= j1")
IF (.NOT. (is_integer(m1) .OR. is_integer(2.0_dp*m1))) &
END IF
IF (.NOT. (is_integer(m1) .OR. is_integer(2.0_dp*m1))) THEN
CPABORT("The angular momentum quantum number m1 has to be integer or half-integer")
IF ((ABS(m2) - j2) > EPSILON(m2)) &
END IF
IF ((ABS(m2) - j2) > EPSILON(m2)) THEN
CPABORT("The angular momentum quantum number m2 has to satisfy -j2 <= m1 <= j2")
IF (.NOT. (is_integer(m2) .OR. is_integer(2.0_dp*m2))) &
END IF
IF (.NOT. (is_integer(m2) .OR. is_integer(2.0_dp*m2))) THEN
CPABORT("The angular momentum quantum number m2 has to be integer or half-integer")
IF ((ABS(M) - J) > EPSILON(M)) &
END IF
IF ((ABS(M) - J) > EPSILON(M)) THEN
CPABORT("The angular momentum quantum number M has to satisfy -J <= M <= J")
IF (.NOT. (is_integer(M) .OR. is_integer(2.0_dp*M))) &
END IF
IF (.NOT. (is_integer(M) .OR. is_integer(2.0_dp*M))) THEN
CPABORT("The angular momentum quantum number M has to be integer or half-integer")
END IF
IF (is_integer(j1 + j2 + J) .AND. &
is_integer(j1 + m1) .AND. &
@ -1872,8 +1884,8 @@ CONTAINS
! **************************************************************************************************
!> \brief Compute the Wigner 3-j symbol
!> / j1 j2 j3 \
!> \ m1 m2 m3 /
!> ( j1 j2 j3 )
!> ( m1 m2 m3 )
!> using the Clebsch-Gordon coefficients
!> \param j1 Angular momentum quantum number of the first state | j1 m1 >
!> \param m1 Magnetic quantum number of the first first state | j1 m1 >

12
src/common/splines.F
View file

@ -311,21 +311,21 @@ CONTAINS
i1 = 1
IF (n < 2) THEN
CALL stop_error("error in iix: n < 2")
ELSEIF (n == 2) THEN
ELSE IF (n == 2) THEN
i1 = 1
ELSEIF (n == 3) THEN
ELSE IF (n == 3) THEN
IF (x <= xi(2)) THEN ! first element
i1 = 1
ELSE
i1 = 2
END IF
ELSEIF (x <= xi(1)) THEN ! left end
ELSE IF (x <= xi(1)) THEN ! left end
i1 = 1
ELSEIF (x <= xi(2)) THEN ! first element
ELSE IF (x <= xi(2)) THEN ! first element
i1 = 1
ELSEIF (x <= xi(3)) THEN ! second element
ELSE IF (x <= xi(3)) THEN ! second element
i1 = 2
ELSEIF (x >= xi(n)) THEN ! right end
ELSE IF (x >= xi(n)) THEN ! right end
i1 = n - 1
ELSE
! bisection: xi(i1) <= x < xi(i2)

21
src/common/timings.F
View file

@ -96,8 +96,9 @@ CONTAINS
IF (PRESENT(timer_env)) timer_env_ => timer_env
IF (.NOT. PRESENT(timer_env)) CALL timer_env_create(timer_env_)
IF (.NOT. ASSOCIATED(timer_env_)) &
IF (.NOT. ASSOCIATED(timer_env_)) THEN
CPABORT("add_timer_env: not associated")
END IF
CALL timer_env_retain(timer_env_)
IF (.NOT. list_isready(timers_stack)) CALL list_init(timers_stack)
@ -157,10 +158,12 @@ CONTAINS
SUBROUTINE timer_env_retain(timer_env)
TYPE(timer_env_type), POINTER :: timer_env
IF (.NOT. ASSOCIATED(timer_env)) &
IF (.NOT. ASSOCIATED(timer_env)) THEN
CPABORT("timer_env_retain: not associated")
IF (timer_env%ref_count < 0) &
END IF
IF (timer_env%ref_count < 0) THEN
CPABORT("timer_env_retain: negativ ref_count")
END IF
timer_env%ref_count = timer_env%ref_count + 1
END SUBROUTINE timer_env_retain
@ -176,10 +179,12 @@ CONTAINS
TYPE(callgraph_item_type), DIMENSION(:), POINTER :: ct_items
TYPE(routine_stat_type), POINTER :: r_stat
IF (.NOT. ASSOCIATED(timer_env)) &
IF (.NOT. ASSOCIATED(timer_env)) THEN
CPABORT("timer_env_release: not associated")
IF (timer_env%ref_count < 0) &
END IF
IF (timer_env%ref_count < 0) THEN
CPABORT("timer_env_release: negativ ref_count")
END IF
timer_env%ref_count = timer_env%ref_count - 1
IF (timer_env%ref_count > 0) RETURN
@ -437,10 +442,12 @@ CONTAINS
TYPE(timer_env_type), POINTER :: timer_env
! catch edge cases where timer_env is not yet/anymore available
IF (.NOT. list_isready(timers_stack)) &
IF (.NOT. list_isready(timers_stack)) THEN
RETURN
IF (list_size(timers_stack) == 0) &
END IF
IF (list_size(timers_stack) == 0) THEN
RETURN
END IF
timer_env => list_peek(timers_stack)
WRITE (unit_nr, '(/,A,/)') " ===== Routine Calling Stack ===== "

9
src/common/timings_report.F
View file

@ -77,8 +77,9 @@ CONTAINS
CALL list_init(reports)
CALL collect_reports_from_ranks(reports, cost_type, para_env)
IF (list_size(reports) > 0 .AND. iw > 0) &
IF (list_size(reports) > 0 .AND. iw > 0) THEN
CALL print_reports(reports, iw, r_timings, sort_by_self_time, cost_type, report_maxloc, para_env)
END IF
! deallocate reports
DO WHILE (list_size(reports) > 0)
@ -111,8 +112,9 @@ CONTAINS
TYPE(timer_env_type), POINTER :: timer_env
NULLIFY (r_stat, r_report, timer_env)
IF (.NOT. list_isready(reports)) &
IF (.NOT. list_isready(reports)) THEN
CPABORT("BUG")
END IF
timer_env => get_timer_env()
@ -232,8 +234,9 @@ CONTAINS
TYPE(routine_report_type), POINTER :: r_report_i, r_report_j
NULLIFY (r_report_i, r_report_j)
IF (.NOT. list_isready(reports)) &
IF (.NOT. list_isready(reports)) THEN
CPABORT("BUG")
END IF
! are we printing timing or energy ?
SELECT CASE (cost_type)

2
src/common/util.F
View file

@ -174,7 +174,7 @@ CONTAINS
IF (j == 1) THEN
DO i = 1, isize
INDEX(i) = i
ENDDO
END DO
END IF
! Allocate scratch arrays

316
src/commutator_rpnl.F
View file

@ -217,7 +217,7 @@ CONTAINS
nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
ELSEIF (spot) THEN
ELSE IF (spot) THEN
CPABORT('SGP not implemented')
ELSE
CPABORT('PPNL unknown')
@ -868,282 +868,304 @@ CONTAINS
IF (my_rxrv) THEN
! x-component (y [z,Vnl] - z [y, Vnl])
! with LAPACK
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 9), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! yzV
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 3), na, &
! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -yVz
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 9), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -zyV
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 4), na, &
! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! zVy
! with MATMUL
IF (iatom <= jatom) THEN
! yzV
blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! -yVz
blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! -yVz
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4)))
! -zyV
blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -zyV
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! zVy
blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! zVy
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 3)))
ELSE
! yzV
blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1)))
! -yVz
blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4))) ! -yVz
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4)))
! -zyV
blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! -zyV
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1)))
! zVy
blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 3))) ! zVy
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 3)))
END IF
! y-component (z [x,Vnl] - x [z, Vnl])
! with LAPACK
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 7), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! zxV
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 4), na, &
! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -zVx
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 7), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -xzV
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 2), na, &
! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! xVz
! with MATMUL
IF (iatom <= jatom) THEN
! zxV
blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zxV
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! -zVx
blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! -zVx
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 2)))
! -xzV
blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -xzV
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! xVz
blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! xVz
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4)))
ELSE
! zxV
blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! zxV
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1)))
! -zVx
blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 2))) ! -zVx
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 2)))
! -xzV
blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! -xzV
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1)))
! xVz
blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4))) ! xVz
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4)))
END IF
! z-component (x [y,Vnl] - y [x, Vnl])
! with LAPACK
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 6), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! xyV
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 2), na, &
! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -xVy
! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 6), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -yxV
! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 3), na, &
! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! yVx
! with MATMUL
IF (iatom <= jatom) THEN
! xyV
blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! -xVy
blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! -xVy
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3)))
! -yxV
blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -yxV
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! zVx
blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! zVx
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 2)))
ELSE
! xyV
blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1)))
! -xVy
blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3))) ! -xVy
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3)))
! -yxV
blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! -yxV
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1)))
! zVx
blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 2))) ! zVx
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 2)))
END IF
END IF
IF (my_rrv) THEN
! r_alpha * r_beta * Vnl
! with LAPACK
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 5), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! xxV
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 6), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! xyV
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 7), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! xzV
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 8), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! yyV
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 9), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! yzV
! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 10), na, &
! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! zzV
! with MATMUL
IF (iatom <= jatom) THEN
! xxV
blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xxV
MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! xyV
blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! xzV
blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xzV
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! yyV
blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yyV
MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! yzV
blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! zzV
blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zzV
MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1)))
ELSE
! xxV
blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1))) ! xxV
MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1)))
! xyV
blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1)))
! xzV
blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! xzV
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1)))
! yyV
blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1))) ! yyV
MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1)))
! yzV
blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1)))
! zzV
blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1))) ! zzV
MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1)))
END IF
! - Vnl * r_alpha * r_beta
! with LAPACK
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 5), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! Vxx
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 6), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! Vxy
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 7), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! Vxz
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 8), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! Vyy
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 9), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! Vyz
! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
! bcint(1, 1, 10), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! Vzz
! with MATMUL
IF (iatom <= jatom) THEN
! -Vxx
blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5))) ! -Vxx
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5)))
! -Vxy
blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6))) ! -Vxy
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6)))
! -Vxz
blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7))) ! -Vxz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7)))
! -Vyy
blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8))) ! -Vyy
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8)))
! -Vyz
blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9))) ! -Vyz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9)))
! -Vzz
blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10))) ! -Vzz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10)))
ELSE
! -Vxx
blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5))) ! -Vxx
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5)))
! -Vxy
blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6))) ! -Vxy
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6)))
! -Vxz
blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7))) ! -Vxz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7)))
! -Vyy
blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8))) ! -Vyy
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8)))
! -Vyz
blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9))) ! -Vyz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9)))
! -Vzz
blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10))) ! -Vzz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10)))
END IF
END IF
IF (my_rvr) THEN
! r_alpha * Vnl * r_beta
IF (iatom <= jatom) THEN
! xVx
blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! xVx
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 2)))
! xVy
blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! xVy
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3)))
! xVz
blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! xVz
MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4)))
! yVy
blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! yVy
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 3)))
! yVz
blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! yVz
MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4)))
! zVz
blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! zVz
MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 4)))
ELSE
! xVx
blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 2))) ! xVx
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 2)))
! xVy
blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3))) ! xVy
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3)))
! xVz
blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4))) ! xVz
MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4)))
! yVy
blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 3))) ! yVy
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 3)))
! yVz
blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4))) ! yVz
MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4)))
! zVz
blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 4))) ! zVz
MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 4)))
END IF
END IF
IF (my_rrv_vrr) THEN
! r_alpha * r_beta * Vnl
IF (iatom <= jatom) THEN
! xxV
blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xxV
MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! xyV
blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! xzV
blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xzV
MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! yyV
blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yyV
MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! yzV
blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1)))
! zzV
blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zzV
MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1)))
ELSE
! xxV
blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1))) ! xxV
MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1)))
! xyV
blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1)))
! xzV
blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! xzV
MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1)))
! yyV
blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1))) ! yyV
MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1)))
! yzV
blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1)))
! zzV
blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1))) ! zzV
MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1)))
END IF
! + Vnl * r_alpha * r_beta
IF (iatom <= jatom) THEN
! +Vxx
blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5))) ! +Vxx
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5)))
! +Vxy
blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6))) ! +Vxy
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6)))
! +Vxz
blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7))) ! +Vxz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7)))
! +Vyy
blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8))) ! +Vyy
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8)))
! +Vyz
blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9))) ! +Vyz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9)))
! +Vzz
blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10))) ! +Vzz
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10)))
ELSE
! +Vxx
blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5))) ! +Vxx
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5)))
! +Vxy
blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6))) ! +Vxy
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6)))
! +Vxz
blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7))) ! +Vxz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7)))
! +Vyy
blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8))) ! +Vyy
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8)))
! +Vyz
blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9))) ! +Vyz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9)))
! +Vzz
blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10))) ! +Vzz
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10)))
END IF
END IF
! The indices are stored in i_1, i_x, ..., i_zzz
! matrix_r_rxvr(alpha, beta)
! = sum_(gamma delta) epsilon_(alpha gamma delta)
! (r_beta * r_gamma * V_nl * r_delta - r_beta * r_gamma * r_delta * V_nl)
! = sum_(gamma delta) epsilon_(alpha gamma delta) r_beta * r_gamma * V_nl * r_delta
! TODO: is this set to zero before?
IF (my_r_rxvr) THEN

95
src/constraint.F
View file

@ -157,43 +157,51 @@ CONTAINS
int_max_sigma = 0.0_dp
ishake_int = ishake_int + 1
! 3x3
IF (n3x3con /= 0) &
IF (n3x3con /= 0) THEN
CALL shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake_int, &
int_max_sigma)
END IF
! 4x6
IF (n4x6con /= 0) &
IF (n4x6con /= 0) THEN
CALL shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake_int, &
int_max_sigma)
END IF
! Collective Variables
IF (ncolv%ntot /= 0) &
IF (ncolv%ntot /= 0) THEN
CALL shake_colv_int(molecule, particle_set, pos, vel, dt, ishake_int, &
cell, imass, int_max_sigma)
END IF
END DO Shake_Intra_Loop
max_sigma = MAX(max_sigma, int_max_sigma)
CALL shake_int_info(log_unit, i, ishake_int, max_sigma)
! Virtual Site
IF (nvsitecon /= 0) &
IF (nvsitecon /= 0) THEN
CALL shake_vsite_int(molecule, pos)
END IF
END DO
END DO MOL
! Intermolecular constraints
IF (do_ext_constraint) THEN
CALL update_temporary_set(group, pos=pos, vel=vel)
! 3x3
IF (gci%ng3x3 /= 0) &
IF (gci%ng3x3 /= 0) THEN
CALL shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
max_sigma)
END IF
! 4x6
IF (gci%ng4x6 /= 0) &
IF (gci%ng4x6 /= 0) THEN
CALL shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
max_sigma)
END IF
! Collective Variables
IF (gci%ncolv%ntot /= 0) &
IF (gci%ncolv%ntot /= 0) THEN
CALL shake_colv_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
cell, imass, max_sigma)
END IF
! Virtual Site
IF (gci%nvsite /= 0) &
IF (gci%nvsite /= 0) THEN
CALL shake_vsite_ext(gci, pos)
END IF
CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel)
END IF
CALL shake_ext_info(log_unit, ishake_ext, max_sigma)
@ -283,15 +291,18 @@ CONTAINS
int_max_sigma = 0.0_dp
irattle_int = irattle_int + 1
! 3x3
IF (n3x3con /= 0) &
IF (n3x3con /= 0) THEN
CALL rattle_3x3_int(molecule, particle_set, vel, dt)
END IF
! 4x6
IF (n4x6con /= 0) &
IF (n4x6con /= 0) THEN
CALL rattle_4x6_int(molecule, particle_set, vel, dt)
END IF
! Collective Variables
IF (ncolv%ntot /= 0) &
IF (ncolv%ntot /= 0) THEN
CALL rattle_colv_int(molecule, particle_set, vel, dt, &
irattle_int, cell, imass, int_max_sigma)
END IF
END DO Rattle_Intra_Loop
max_sigma = MAX(max_sigma, int_max_sigma)
CALL rattle_int_info(log_unit, i, irattle_int, max_sigma)
@ -301,15 +312,18 @@ CONTAINS
IF (do_ext_constraint) THEN
CALL update_temporary_set(group, vel=vel)
! 3x3
IF (gci%ng3x3 /= 0) &
IF (gci%ng3x3 /= 0) THEN
CALL rattle_3x3_ext(gci, particle_set, vel, dt)
END IF
! 4x6
IF (gci%ng4x6 /= 0) &
IF (gci%ng4x6 /= 0) THEN
CALL rattle_4x6_ext(gci, particle_set, vel, dt)
END IF
! Collective Variables
IF (gci%ncolv%ntot /= 0) &
IF (gci%ncolv%ntot /= 0) THEN
CALL rattle_colv_ext(gci, particle_set, vel, dt, &
irattle_ext, cell, imass, max_sigma)
END IF
CALL restore_temporary_set(particle_set, local_particles, vel=vel)
END IF
CALL rattle_ext_info(log_unit, irattle_ext, max_sigma)
@ -417,17 +431,20 @@ CONTAINS
int_max_sigma = 0.0_dp
ishake_int = ishake_int + 1
! 3x3
IF (n3x3con /= 0) &
IF (n3x3con /= 0) THEN
CALL shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, &
v_shake, dt, ishake_int, int_max_sigma)
END IF
! 4x6
IF (n4x6con /= 0) &
IF (n4x6con /= 0) THEN
CALL shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, &
dt, ishake_int, int_max_sigma)
END IF
! Collective Variables
IF (ncolv%ntot /= 0) &
IF (ncolv%ntot /= 0) THEN
CALL shake_roll_colv_int(molecule, particle_set, pos, vel, r_shake, &
v_shake, dt, ishake_int, cell, imass, int_max_sigma)
END IF
END DO Shake_Roll_Intra_Loop
max_sigma = MAX(max_sigma, int_max_sigma)
CALL shake_int_info(log_unit, i, ishake_int, max_sigma)
@ -441,20 +458,24 @@ CONTAINS
IF (do_ext_constraint) THEN
CALL update_temporary_set(group, pos=pos, vel=vel)
! 3x3
IF (gci%ng3x3 /= 0) &
IF (gci%ng3x3 /= 0) THEN
CALL shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, &
v_shake, dt, ishake_ext, max_sigma)
END IF
! 4x6
IF (gci%ng4x6 /= 0) &
IF (gci%ng4x6 /= 0) THEN
CALL shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, &
dt, ishake_ext, max_sigma)
END IF
! Collective Variables
IF (gci%ncolv%ntot /= 0) &
IF (gci%ncolv%ntot /= 0) THEN
CALL shake_roll_colv_ext(gci, particle_set, pos, vel, r_shake, &
v_shake, dt, ishake_ext, cell, imass, max_sigma)
END IF
! Virtual Site
IF (gci%nvsite /= 0) &
IF (gci%nvsite /= 0) THEN
CPABORT("Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!")
END IF
CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel)
END IF
CALL shake_ext_info(log_unit, ishake_ext, max_sigma)
@ -564,17 +585,20 @@ CONTAINS
int_max_sigma = 0.0_dp
irattle_int = irattle_int + 1
! 3x3
IF (n3x3con /= 0) &
IF (n3x3con /= 0) THEN
CALL rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, &
veps)
END IF
! 4x6
IF (n4x6con /= 0) &
IF (n4x6con /= 0) THEN
CALL rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, &
veps)
END IF
! Collective Variables
IF (ncolv%ntot /= 0) &
IF (ncolv%ntot /= 0) THEN
CALL rattle_roll_colv_int(molecule, particle_set, vel, r_rattle, dt, &
irattle_int, veps, cell, imass, int_max_sigma)
END IF
END DO Rattle_Roll_Intramolecular
max_sigma = MAX(max_sigma, int_max_sigma)
CALL rattle_int_info(log_unit, i, irattle_int, max_sigma)
@ -584,17 +608,20 @@ CONTAINS
IF (do_ext_constraint) THEN
CALL update_temporary_set(para_env, vel=vel)
! 3x3
IF (gci%ng3x3 /= 0) &
IF (gci%ng3x3 /= 0) THEN
CALL rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, &
veps)
END IF
! 4x6
IF (gci%ng4x6 /= 0) &
IF (gci%ng4x6 /= 0) THEN
CALL rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, &
veps)
END IF
! Collective Variables
IF (gci%ncolv%ntot /= 0) &
IF (gci%ncolv%ntot /= 0) THEN
CALL rattle_roll_colv_ext(gci, particle_set, vel, r_rattle, dt, &
irattle_ext, veps, cell, imass, max_sigma)
END IF
CALL restore_temporary_set(particle_set, local_particles, vel=vel)
END IF
CALL rattle_ext_info(log_unit, irattle_ext, max_sigma)
@ -708,7 +735,7 @@ CONTAINS
IF (log_unit > 0) THEN
IF (id_type == "S") THEN
label = "Shake Lagrangian Multipliers:"
ELSEIF (id_type == "R") THEN
ELSE IF (id_type == "R") THEN
label = "Rattle Lagrangian Multipliers:"
ELSE
CPABORT("Only S for Shake or R for Rattle are supported for Lagrangian Multipliers")
@ -743,11 +770,12 @@ CONTAINS
"Molecule Nr.:", i, " Nr. Iterations:", ishake_int, " Max. Err.:", max_sigma
END IF
! Notify a not converged SHAKE
IF (ishake_int > Max_Shake_Iter) &
IF (ishake_int > Max_Shake_Iter) THEN
CALL cp_warn(__LOCATION__, &
"Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
"intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// &
". CP2K continues but results could be meaningless. ")
END IF
END SUBROUTINE shake_int_info
! **************************************************************************************************
@ -769,10 +797,11 @@ CONTAINS
" Max. Err.:", max_sigma
END IF
! Notify a not converged SHAKE
IF (ishake_ext > Max_Shake_Iter) &
IF (ishake_ext > Max_Shake_Iter) THEN
CALL cp_warn(__LOCATION__, &
"Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
"intermolecular constraint. CP2K continues but results could be meaningless.")
END IF
END SUBROUTINE shake_ext_info
! **************************************************************************************************
@ -794,11 +823,12 @@ CONTAINS
"Molecule Nr.:", i, " Nr. Iterations:", irattle_int, " Max. Err.:", max_sigma
END IF
! Notify a not converged RATTLE
IF (irattle_int > Max_shake_Iter) &
IF (irattle_int > Max_shake_Iter) THEN
CALL cp_warn(__LOCATION__, &
"Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
"intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// &
". CP2K continues but results could be meaningless.")
END IF
END SUBROUTINE rattle_int_info
! **************************************************************************************************
@ -820,10 +850,11 @@ CONTAINS
" Max. Err.:", max_sigma
END IF
! Notify a not converged RATTLE
IF (irattle_ext > Max_shake_Iter) &
IF (irattle_ext > Max_shake_Iter) THEN
CALL cp_warn(__LOCATION__, &
"Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
"intermolecular constraint. CP2K continues but results could be meaningless.")
END IF
END SUBROUTINE rattle_ext_info
! **************************************************************************************************

3
src/constraint_fxd.F
View file

@ -398,8 +398,9 @@ CONTAINS
DO k = 1, SIZE(fixd_list)
IF (fixd_list(k)%fixd == j) THEN
IF (fixd_list(k)%itype /= use_perd_xyz) CYCLE
IF (.NOT. fixd_list(k)%restraint%active) &
IF (.NOT. fixd_list(k)%restraint%active) THEN
colvar%dsdr(:, i) = 0.0_dp
END IF
EXIT
END IF
END DO

6
src/constraint_util.F
View file

@ -554,8 +554,8 @@ CONTAINS
v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
diag = MATMUL(r_shake, v_shake)
r_shake = diag
ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. &
PRESENT(vector_r)) THEN
ELSE IF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. &
PRESENT(vector_r)) THEN
DO i = 1, 3
r_shake(i, i) = vector_r(i)*vector_v(i)
v_shake(i, i) = vector_v(i)
@ -569,7 +569,7 @@ CONTAINS
diag(2, 2) = vector_v(2)
diag(3, 3) = vector_v(3)
v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN
ELSE IF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN
DO i = 1, 3
v_shake(i, i) = vector_v(i)
END DO

6
src/constraint_vsite.F
View file

@ -94,14 +94,16 @@ CONTAINS
molecule_kind => molecule%molecule_kind
CALL get_molecule_kind(molecule_kind, nconstraint=nconstraint, nvsite=nvsitecon)
IF (nconstraint == 0) CYCLE
IF (nvsitecon /= 0) &
IF (nvsitecon /= 0) THEN
CALL force_vsite_int(molecule, particle_set)
END IF
END DO
END DO MOL
! Intermolecular Virtual Site Constraints
IF (do_ext_constraint) THEN
IF (gci%nvsite /= 0) &
IF (gci%nvsite /= 0) THEN
CALL force_vsite_ext(gci, particle_set)
END IF
END IF
END SUBROUTINE vsite_force_control

8
src/core_ppl.F
View file

@ -448,7 +448,7 @@ CONTAINS
IF (ecp_semi_local) THEN
CALL get_potential(potential=sgp_potential, sl_lmax=slmax, &
npot=npot, nrpot=nrpot, apot=apot, bpot=bpot)
ELSEIF (ecp_local) THEN
ELSE IF (ecp_local) THEN
IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
END IF
ELSE
@ -498,7 +498,7 @@ CONTAINS
rab, dab, rac, dac, rbc, dbc, &
hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
force_a, force_b, ppl_fwork)
ELSEIF (libgrpp_local) THEN
ELSE IF (libgrpp_local) THEN
!$OMP CRITICAL(type1)
CALL libgrpp_local_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
rpgfa(:, iset), zeta(:, iset), &
@ -556,7 +556,7 @@ CONTAINS
CALL virial_pair_force(pv_thread, f0, force_a, rac)
CALL virial_pair_force(pv_thread, f0, force_b, rbc)
END IF
ELSEIF (do_dR) THEN
ELSE IF (do_dR) THEN
hab2_w = 0._dp
CALL ppl_integral( &
la_max(iset), la_min(iset), npgfa(iset), &
@ -584,7 +584,7 @@ CONTAINS
nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
ELSEIF (libgrpp_local) THEN
ELSE IF (libgrpp_local) THEN
!If the local part of the potential is more complex, we need libgrpp
!$OMP CRITICAL(type1)
CALL libgrpp_local_integrals(la_max(iset), la_min(iset), npgfa(iset), &

6
src/core_ppnl.F
View file

@ -667,7 +667,7 @@ CONTAINS
END IF
IF (do_dR) THEN
i = 1; j = 2;
i = 1; j = 2
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
@ -686,7 +686,7 @@ CONTAINS
MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
END IF
i = 2; j = 3;
i = 2; j = 3
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
@ -705,7 +705,7 @@ CONTAINS
MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
END IF
i = 3; j = 4;
i = 3; j = 4
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &

2
src/cp2k_debug.F
View file

@ -612,7 +612,7 @@ CONTAINS
IF (dft_control%apply_efield) THEN
dft_control%efield_fields(1)%efield%strength = amplitude
dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
ELSEIF (dft_control%apply_period_efield) THEN
ELSE IF (dft_control%apply_period_efield) THEN
dft_control%period_efield%strength = amplitude
dft_control%period_efield%polarisation(1:3) = poldir(1:3)
ELSE

21
src/cp_control_types.F
View file

@ -891,8 +891,9 @@ CONTAINS
SUBROUTINE mulliken_control_release(mulliken_restraint_control)
TYPE(mulliken_restraint_type), INTENT(INOUT) :: mulliken_restraint_control
IF (ASSOCIATED(mulliken_restraint_control%atoms)) &
IF (ASSOCIATED(mulliken_restraint_control%atoms)) THEN
DEALLOCATE (mulliken_restraint_control%atoms)
END IF
mulliken_restraint_control%strength = 0.0_dp
mulliken_restraint_control%target = 0.0_dp
mulliken_restraint_control%natoms = 0
@ -927,10 +928,12 @@ CONTAINS
SUBROUTINE ddapc_control_release(ddapc_restraint_control)
TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
IF (ASSOCIATED(ddapc_restraint_control%atoms)) &
IF (ASSOCIATED(ddapc_restraint_control%atoms)) THEN
DEALLOCATE (ddapc_restraint_control%atoms)
IF (ASSOCIATED(ddapc_restraint_control%coeff)) &
END IF
IF (ASSOCIATED(ddapc_restraint_control%coeff)) THEN
DEALLOCATE (ddapc_restraint_control%coeff)
END IF
ddapc_restraint_control%strength = 0.0_dp
ddapc_restraint_control%target = 0.0_dp
ddapc_restraint_control%natoms = 0
@ -1207,18 +1210,21 @@ CONTAINS
IF (ASSOCIATED(proj_mo_list)) THEN
DO i = 1, SIZE(proj_mo_list)
IF (ASSOCIATED(proj_mo_list(i)%proj_mo)) THEN
IF (ALLOCATED(proj_mo_list(i)%proj_mo%ref_mo_index)) &
IF (ALLOCATED(proj_mo_list(i)%proj_mo%ref_mo_index)) THEN
DEALLOCATE (proj_mo_list(i)%proj_mo%ref_mo_index)
END IF
IF (ALLOCATED(proj_mo_list(i)%proj_mo%mo_ref)) THEN
DO mo_ref_nbr = 1, SIZE(proj_mo_list(i)%proj_mo%mo_ref)
CALL cp_fm_release(proj_mo_list(i)%proj_mo%mo_ref(mo_ref_nbr))
END DO
DEALLOCATE (proj_mo_list(i)%proj_mo%mo_ref)
END IF
IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_index)) &
IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_index)) THEN
DEALLOCATE (proj_mo_list(i)%proj_mo%td_mo_index)
IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_occ)) &
END IF
IF (ALLOCATED(proj_mo_list(i)%proj_mo%td_mo_occ)) THEN
DEALLOCATE (proj_mo_list(i)%proj_mo%td_mo_occ)
END IF
DEALLOCATE (proj_mo_list(i)%proj_mo)
END IF
END DO
@ -1244,8 +1250,9 @@ CONTAINS
IF (ASSOCIATED(efield_fields(i)%efield%envelop_i_vars)) THEN
DEALLOCATE (efield_fields(i)%efield%envelop_i_vars)
END IF
IF (ASSOCIATED(efield_fields(i)%efield%polarisation)) &
IF (ASSOCIATED(efield_fields(i)%efield%polarisation)) THEN
DEALLOCATE (efield_fields(i)%efield%polarisation)
END IF
DEALLOCATE (efield_fields(i)%efield)
END IF
END DO

149
src/cp_control_utils.F
View file

@ -164,20 +164,23 @@ CONTAINS
CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut)
CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut)
! Perform numerical stability checks and possibly correct the issues
IF (density_cut <= EPSILON(0.0_dp)*100.0_dp) &
IF (density_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
CALL cp_warn(__LOCATION__, &
"DENSITY_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
"This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
END IF
density_cut = MAX(EPSILON(0.0_dp)*100.0_dp, density_cut)
IF (gradient_cut <= EPSILON(0.0_dp)*100.0_dp) &
IF (gradient_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
CALL cp_warn(__LOCATION__, &
"GRADIENT_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
"This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
END IF
gradient_cut = MAX(EPSILON(0.0_dp)*100.0_dp, gradient_cut)
IF (tau_cut <= EPSILON(0.0_dp)*100.0_dp) &
IF (tau_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
CALL cp_warn(__LOCATION__, &
"TAU_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
"This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
END IF
tau_cut = MAX(EPSILON(0.0_dp)*100.0_dp, tau_cut)
CALL section_vals_val_set(xc_section, "density_cutoff", r_val=density_cut)
CALL section_vals_val_set(xc_section, "gradient_cutoff", r_val=gradient_cut)
@ -416,8 +419,9 @@ CONTAINS
IF (dft_control%admm_control%purification_method == do_admm_purify_mo_diag .OR. &
dft_control%admm_control%purification_method == do_admm_purify_mo_no_diag) THEN
IF (dft_control%admm_control%method /= do_admm_basis_projection) &
IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
CPABORT("ADMM: Chosen purification requires BASIS_PROJECTION")
END IF
IF (.NOT. do_ot) CPABORT("ADMM: MO-based purification requires OT.")
END IF
@ -440,9 +444,10 @@ CONTAINS
CALL section_vals_val_get(dft_section, "MULTIPLICITY", i_val=dft_control%multiplicity)
CALL section_vals_val_get(dft_section, "RELAX_MULTIPLICITY", r_val=dft_control%relax_multiplicity)
IF (dft_control%relax_multiplicity > 0.0_dp) THEN
IF (.NOT. dft_control%uks) &
IF (.NOT. dft_control%uks) THEN
CALL cp_abort(__LOCATION__, "The option RELAX_MULTIPLICITY is only valid for "// &
"unrestricted Kohn-Sham (UKS) calculations")
END IF
END IF
!Read the HAIR PROBES input section if present
@ -578,8 +583,9 @@ CONTAINS
CALL section_vals_val_get(tmp_section, "POLARISATION", r_vals=pol)
dft_control%period_efield%polarisation(1:3) = pol(1:3)
IF (PRESENT(cell)) THEN
IF (ASSOCIATED(cell)) &
IF (ASSOCIATED(cell)) THEN
CALL cell_transform_input_cartesian(cell, dft_control%period_efield%polarisation(1:3))
END IF
END IF
CALL section_vals_val_get(tmp_section, "D_FILTER", r_vals=pol)
dft_control%period_efield%d_filter(1:3) = pol(1:3)
@ -659,11 +665,12 @@ CONTAINS
END IF
! periodic fields don't work with RTP
IF (do_rtp) &
IF (do_rtp) THEN
CALL cp_abort(__LOCATION__, &
"Periodic efield cannot be used with RTP. When restarting a "// &
"run with periodic efield, set RESTART_RTP under &EXT_RESTART "// &
"section to .FALSE. explicitly if RESTART_DEFAULT is .TRUE.")
END IF
IF (dft_control%period_efield%displacement_field) THEN
CALL cite_reference(Stengel2009)
ELSE
@ -1232,8 +1239,9 @@ CONTAINS
jj = jj + SIZE(tmplist)
END DO
qs_control%mulliken_restraint_control%natoms = jj
IF (qs_control%mulliken_restraint_control%natoms < 1) &
IF (qs_control%mulliken_restraint_control%natoms < 1) THEN
CPABORT("Need at least 1 atom to use mulliken constraints")
END IF
ALLOCATE (qs_control%mulliken_restraint_control%atoms(qs_control%mulliken_restraint_control%natoms))
jj = 0
DO k = 1, n_rep
@ -1281,10 +1289,11 @@ CONTAINS
CALL section_vals_val_get(se_section, "INTEGRAL_SCREENING", &
i_val=qs_control%se_control%integral_screening)
IF (qs_control%method_id == do_method_pnnl) THEN
IF (qs_control%se_control%integral_screening /= do_se_IS_slater) &
IF (qs_control%se_control%integral_screening /= do_se_IS_slater) THEN
CALL cp_warn(__LOCATION__, &
"PNNL semi-empirical parameterization supports only the Slater type "// &
"integral scheme. Revert to Slater and continue the calculation.")
END IF
qs_control%se_control%integral_screening = do_se_IS_slater
END IF
! Global Arrays variable
@ -1349,21 +1358,23 @@ CONTAINS
qs_control%se_control%do_ewald = .FALSE.
qs_control%se_control%do_ewald_r3 = .FALSE.
qs_control%se_control%do_ewald_gks = .TRUE.
IF (qs_control%method_id /= do_method_pnnl) &
IF (qs_control%method_id /= do_method_pnnl) THEN
CALL cp_abort(__LOCATION__, &
"A periodic semi-empirical calculation was requested with a long-range "// &
"summation on the single integral evaluation. This scheme is supported "// &
"only by the PNNL parameterization.")
END IF
CASE (do_se_lr_ewald_r3)
qs_control%se_control%do_ewald = .TRUE.
qs_control%se_control%do_ewald_r3 = .TRUE.
qs_control%se_control%do_ewald_gks = .FALSE.
IF (qs_control%se_control%integral_screening /= do_se_IS_kdso) &
IF (qs_control%se_control%integral_screening /= do_se_IS_kdso) THEN
CALL cp_abort(__LOCATION__, &
"A periodic semi-empirical calculation was requested with a long-range "// &
"summation for the slowly convergent part 1/R^3, which is not congruent "// &
"with the integral screening chosen. The only integral screening supported "// &
"by this periodic type calculation is the standard Klopman-Dewar-Sabelli-Ohno.")
END IF
END SELECT
! dispersion pair potentials
@ -1427,18 +1438,21 @@ CONTAINS
"DFTB/TBLITE_MIXER")
IF (qs_control%do_ls_scf) THEN
IF (dftb_scc_mixer_explicit .AND. &
qs_control%dftb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
qs_control%dftb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
CALL cp_warn(__LOCATION__, &
"DFTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
"the density matrix directly.")
IF (dftb_tblite_mixer_explicit) &
END IF
IF (dftb_tblite_mixer_explicit) THEN
CALL cp_warn(__LOCATION__, &
"DFTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
"the density-matrix optimization.")
END IF
qs_control%dftb_control%tblite_scc_mixer = tblite_scc_mixer_none
END IF
IF (qs_control%dftb_control%tblite_mixer_damping <= 0.0_dp) &
IF (qs_control%dftb_control%tblite_mixer_damping <= 0.0_dp) THEN
CPABORT("DFTB/TBLITE_MIXER/DAMPING must be positive")
END IF
CALL section_vals_val_get(dftb_section, "EPS_DISP", &
r_val=qs_control%dftb_control%eps_disp)
CALL section_vals_val_get(dftb_section, "DO_EWALD", explicit=explicit)
@ -1498,8 +1512,9 @@ CONTAINS
CALL section_vals_val_get(xtb_tblite, "_SECTION_PARAMETERS_", l_val=tblite_section_active)
qs_control%xtb_control%do_tblite = (qs_control%xtb_control%gfn_type == gfn_tblite)
IF (qs_control%xtb_control%do_tblite) THEN
IF (.NOT. tblite_section_active) &
IF (.NOT. tblite_section_active) THEN
CPABORT("XTB/GFN_TYPE TBLITE requires an XTB/TBLITE section")
END IF
! The CP2K-internal GFN1 defaults are still used to initialize shared xTB fields.
qs_control%xtb_control%gfn_type = gfn1xtb
ELSE IF (tblite_section_active) THEN
@ -1520,37 +1535,43 @@ CONTAINS
qs_control%xtb_control%tblite_mixer_max_weight, &
qs_control%xtb_control%tblite_mixer_weight_factor, &
"XTB/TBLITE_MIXER")
IF (xtb_tblite_mixer_explicit) &
IF (xtb_tblite_mixer_explicit) THEN
CALL section_vals_val_get(xtb_tblite_mixer, "DAMPING", &
explicit=qs_control%xtb_control%tblite_mixer_damping_explicit)
END IF
IF ((.NOT. qs_control%xtb_control%do_tblite) .AND. &
qs_control%xtb_control%gfn_type == 0) THEN
IF (xtb_scc_mixer_explicit .AND. &
qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_auto .AND. &
qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
CALL cp_warn(__LOCATION__, &
"XTB/SCC_MIXER is reset to NONE for CP2K-internal GFN0-xTB; "// &
"GFN0-xTB has no SCC variables to mix.")
IF (xtb_tblite_mixer_explicit) &
END IF
IF (xtb_tblite_mixer_explicit) THEN
CALL cp_warn(__LOCATION__, &
"XTB/TBLITE_MIXER settings are ignored for CP2K-internal GFN0-xTB; "// &
"GFN0-xTB has no SCC variables to mix.")
END IF
qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
END IF
IF (qs_control%do_ls_scf) THEN
IF (xtb_scc_mixer_explicit .AND. &
qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
CALL cp_warn(__LOCATION__, &
"XTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
"the density matrix directly.")
IF (xtb_tblite_mixer_explicit) &
END IF
IF (xtb_tblite_mixer_explicit) THEN
CALL cp_warn(__LOCATION__, &
"XTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
"the density-matrix optimization.")
END IF
qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
END IF
IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) &
IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) THEN
CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive")
END IF
CALL section_vals_val_get(xtb_section, "DO_EWALD", explicit=explicit)
IF (explicit) THEN
CALL section_vals_val_get(xtb_section, "DO_EWALD", &
@ -1870,14 +1891,17 @@ CONTAINS
c_val=qs_control%xtb_control%tblite_param_file)
CALL section_vals_val_get(xtb_tblite, "ACCURACY", &
r_val=qs_control%xtb_control%tblite_accuracy)
IF (qs_control%xtb_control%tblite_accuracy <= 0.0_dp) &
IF (qs_control%xtb_control%tblite_accuracy <= 0.0_dp) THEN
CPABORT("XTB/TBLITE/ACCURACY must be positive")
IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) &
END IF
IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) THEN
CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive")
END IF
CALL section_vals_val_get(xtb_tblite, "REFERENCE_CLI", l_val=tblite_reference_cli)
CALL section_vals_get(xtb_tblite_ref_cli, explicit=tblite_reference_cli_section)
IF (tblite_reference_cli .AND. (.NOT. tblite_reference_cli_section)) &
IF (tblite_reference_cli .AND. (.NOT. tblite_reference_cli_section)) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI keyword requires an XTB/TBLITE/REFERENCE_CLI section")
END IF
IF (tblite_reference_cli .OR. tblite_reference_cli_section) THEN
CALL read_xtb_reference_cli_section(xtb_tblite_ref_cli, qs_control%xtb_control%reference_cli, cell)
qs_control%xtb_control%reference_cli%enabled = .TRUE.
@ -1941,8 +1965,9 @@ CONTAINS
IF (omega0 <= 0.0_dp) CPABORT(TRIM(section_name)//"/OMEGA0 must be positive")
IF (min_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MIN_WEIGHT must be positive")
IF (max_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MAX_WEIGHT must be positive")
IF (max_weight < min_weight) &
IF (max_weight < min_weight) THEN
CPABORT(TRIM(section_name)//"/MAX_WEIGHT must not be smaller than MIN_WEIGHT")
END IF
IF (weight_factor <= 0.0_dp) CPABORT(TRIM(section_name)//"/WEIGHT_FACTOR must be positive")
END SUBROUTINE read_tblite_mixer_section
@ -1991,11 +2016,13 @@ CONTAINS
CALL section_vals_val_get(solvation_section, "SOLVENT", c_val=ref_cli%solvation_solvent)
CALL section_vals_val_get(solvation_section, "BORN_KERNEL", i_val=ref_cli%solvation_born_kernel)
CALL section_vals_val_get(solvation_section, "SOLUTION_STATE", i_val=ref_cli%solvation_state)
IF (LEN_TRIM(ref_cli%solvation_solvent) == 0) &
IF (LEN_TRIM(ref_cli%solvation_solvent) == 0) THEN
CPABORT("REFERENCE_CLI implicit solvation needs SOLVENT")
END IF
IF (ref_cli%solvation_model == tblite_cli_solvation_cpcm .AND. &
ref_cli%solvation_born_kernel /= tblite_cli_born_kernel_auto) &
ref_cli%solvation_born_kernel /= tblite_cli_born_kernel_auto) THEN
CPABORT("BORN_KERNEL is invalid with MODEL CPCM")
END IF
IF (ref_cli%solvation_state /= tblite_cli_solution_state_gsolv) THEN
SELECT CASE (ref_cli%solvation_model)
CASE (tblite_cli_solvation_alpb, tblite_cli_solvation_gbsa)
@ -2007,18 +2034,21 @@ CONTAINS
END IF
CALL section_vals_val_get(ref_cli_section, "ELECTRONIC_TEMPERATURE_GUESS", &
r_val=ref_cli%electronic_temperature_guess)
IF (ref_cli%electronic_temperature_guess < 0.0_dp) &
IF (ref_cli%electronic_temperature_guess < 0.0_dp) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
IF (ref_cli%electronic_temperature_guess > 0.0_dp .AND. ref_cli%guess /= tblite_guess_ceh) &
END IF
IF (ref_cli%electronic_temperature_guess > 0.0_dp .AND. ref_cli%guess /= tblite_guess_ceh) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS requires GUESS CEH")
END IF
guess_section => section_vals_get_subs_vals(ref_cli_section, "GUESS_CLI")
CALL section_vals_get(guess_section, explicit=ref_cli%guess_cli%enabled)
IF (ref_cli%guess_cli%enabled) THEN
CALL section_vals_val_get(guess_section, "METHOD", i_val=ref_cli%guess_cli%method)
CALL section_vals_val_get(guess_section, "ELECTRONIC_TEMPERATURE_GUESS", &
r_val=ref_cli%guess_cli%electronic_temperature_guess)
IF (ref_cli%guess_cli%electronic_temperature_guess < 0.0_dp) &
IF (ref_cli%guess_cli%electronic_temperature_guess < 0.0_dp) THEN
CPABORT("REFERENCE_CLI/GUESS_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
END IF
CALL section_vals_val_get(guess_section, "SOLVER", i_val=ref_cli%guess_cli%solver)
CALL section_vals_val_get(guess_section, "EFIELD", explicit=ref_cli%guess_cli%efield_active)
IF (ref_cli%guess_cli%efield_active) THEN
@ -2049,10 +2079,12 @@ CONTAINS
CALL section_vals_val_get(fit_section, "INPUT_FILE", c_val=ref_cli%fit_cli%input_file)
CALL section_vals_val_get(fit_section, "DRY_RUN", l_val=ref_cli%fit_cli%dry_run)
CALL section_vals_val_get(fit_section, "COPY", c_val=ref_cli%fit_cli%copy_file)
IF (LEN_TRIM(ref_cli%fit_cli%param_file) == 0) &
IF (LEN_TRIM(ref_cli%fit_cli%param_file) == 0) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs PARAM_FILE")
IF (LEN_TRIM(ref_cli%fit_cli%input_file) == 0) &
END IF
IF (LEN_TRIM(ref_cli%fit_cli%input_file) == 0) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs INPUT_FILE")
END IF
END IF
tagdiff_section => section_vals_get_subs_vals(ref_cli_section, "TAGDIFF_CLI")
CALL section_vals_get(tagdiff_section, explicit=ref_cli%tagdiff_cli%enabled)
@ -2060,10 +2092,12 @@ CONTAINS
CALL section_vals_val_get(tagdiff_section, "ACTUAL", c_val=ref_cli%tagdiff_cli%actual_file)
CALL section_vals_val_get(tagdiff_section, "REFERENCE", c_val=ref_cli%tagdiff_cli%reference_file)
CALL section_vals_val_get(tagdiff_section, "FIT", l_val=ref_cli%tagdiff_cli%fit)
IF (LEN_TRIM(ref_cli%tagdiff_cli%actual_file) == 0) &
IF (LEN_TRIM(ref_cli%tagdiff_cli%actual_file) == 0) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs ACTUAL")
IF (LEN_TRIM(ref_cli%tagdiff_cli%reference_file) == 0) &
END IF
IF (LEN_TRIM(ref_cli%tagdiff_cli%reference_file) == 0) THEN
CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs REFERENCE")
END IF
END IF
CALL section_vals_val_get(ref_cli_section, "KEEP_FILES", l_val=ref_cli%keep_files)
CALL section_vals_val_get(ref_cli_section, "ERROR_LIMIT", r_val=ref_cli%error_limit)
@ -2174,8 +2208,9 @@ CONTAINS
END IF
END DO
IF (t_control%conv < 0) &
IF (t_control%conv < 0) THEN
t_control%conv = ABS(t_control%conv)
END IF
! DIPOLE_MOMENTS subsection
dipole_section => section_vals_get_subs_vals(t_section, "DIPOLE_MOMENTS")
@ -2217,9 +2252,10 @@ CONTAINS
CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", &
r_val=t_control%mgrid_progression_factor, explicit=explicit)
IF (explicit) THEN
IF (t_control%mgrid_progression_factor <= 1.0_dp) &
IF (t_control%mgrid_progression_factor <= 1.0_dp) THEN
CALL cp_abort(__LOCATION__, &
"Progression factor should be greater then 1.0 to ensure multi-grid ordering")
END IF
ELSE
t_control%mgrid_progression_factor = qs_control%progression_factor
END IF
@ -2253,8 +2289,9 @@ CONTAINS
IF (.NOT. explicit) t_control%mgrid_skip_load_balance = qs_control%skip_load_balance_distributed
IF (ASSOCIATED(t_control%mgrid_e_cutoff)) THEN
IF (SIZE(t_control%mgrid_e_cutoff) /= t_control%mgrid_ngrids) &
IF (SIZE(t_control%mgrid_e_cutoff) /= t_control%mgrid_ngrids) THEN
CPABORT("Inconsistent values for number of multi-grids")
END IF
! sort multi-grids in descending order according to their cutoff values
t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff
@ -2269,8 +2306,9 @@ CONTAINS
xc_section => section_vals_get_subs_vals(t_section, "XC")
xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
CALL section_vals_get(xc_func, explicit=explicit)
IF (explicit) &
IF (explicit) THEN
CALL xc_functionals_expand(xc_func, xc_section)
END IF
! sTDA subsection
stda_section => section_vals_get_subs_vals(t_section, "STDA")
@ -2490,8 +2528,7 @@ CONTAINS
dft_control%period_efield%strength
END IF
IF (SQRT(DOT_PRODUCT(dft_control%period_efield%polarisation, &
dft_control%period_efield%polarisation)) < EPSILON(0.0_dp)) THEN
IF (NORM2(dft_control%period_efield%polarisation) < EPSILON(0.0_dp)) THEN
CPABORT("Invalid (too small) polarisation vector specified for PERIODIC_EFIELD")
END IF
END IF
@ -2900,8 +2937,9 @@ CONTAINS
ELSE
WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
"QS| Density cutoff [a.u.]:", qs_control%cutoff
IF (qs_control%commensurate_mgrids) &
IF (qs_control%commensurate_mgrids) THEN
WRITE (UNIT=output_unit, FMT="(T2,A)") "QS| Using commensurate multigrids"
END IF
WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
"QS| Multi grid cutoff [a.u.]: 1) grid level", qs_control%e_cutoff(1)
WRITE (UNIT=output_unit, FMT="(T2,A,I3,A,T71,F10.1)") &
@ -3031,9 +3069,10 @@ CONTAINS
IF (qs_control%ddapc_restraint) THEN
DO i = 1, SIZE(qs_control%ddapc_restraint_control)
ddapc_restraint_control => qs_control%ddapc_restraint_control(i)
IF (SIZE(qs_control%ddapc_restraint_control) > 1) &
IF (SIZE(qs_control%ddapc_restraint_control) > 1) THEN
WRITE (UNIT=output_unit, FMT="(T2,A,T3,I8)") &
"QS| parameters for DDAPC restraint number", i
"QS| parameters for DDAPC restraint number", i
END IF
WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
"QS| ddapc restraint target", ddapc_restraint_control%target
WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
@ -3104,8 +3143,9 @@ CONTAINS
IF (PRESENT(ddapc_restraint_section)) THEN
IF (ASSOCIATED(qs_control%ddapc_restraint_control)) THEN
IF (SIZE(qs_control%ddapc_restraint_control) >= 2) &
IF (SIZE(qs_control%ddapc_restraint_control) >= 2) THEN
CPABORT("ET_COUPLING cannot be used in combination with a normal restraint")
END IF
ELSE
ddapc_section => ddapc_restraint_section
ALLOCATE (qs_control%ddapc_restraint_control(1))
@ -3144,8 +3184,9 @@ CONTAINS
END DO
IF (jj < 1) CPABORT("Need at least 1 atom to use ddapc constraints")
ddapc_restraint_control%natoms = jj
IF (ASSOCIATED(ddapc_restraint_control%atoms)) &
IF (ASSOCIATED(ddapc_restraint_control%atoms)) THEN
DEALLOCATE (ddapc_restraint_control%atoms)
END IF
ALLOCATE (ddapc_restraint_control%atoms(ddapc_restraint_control%natoms))
jj = 0
DO k = 1, n_rep
@ -3157,8 +3198,9 @@ CONTAINS
END DO
END DO
IF (ASSOCIATED(ddapc_restraint_control%coeff)) &
IF (ASSOCIATED(ddapc_restraint_control%coeff)) THEN
DEALLOCATE (ddapc_restraint_control%coeff)
END IF
ALLOCATE (ddapc_restraint_control%coeff(ddapc_restraint_control%natoms))
ddapc_restraint_control%coeff = 1.0_dp
@ -3170,13 +3212,15 @@ CONTAINS
i_rep_val=k, r_vals=rtmplist)
DO j = 1, SIZE(rtmplist)
jj = jj + 1
IF (jj > ddapc_restraint_control%natoms) &
IF (jj > ddapc_restraint_control%natoms) THEN
CPABORT("Need the same number of coeff as there are atoms ")
END IF
ddapc_restraint_control%coeff(jj) = rtmplist(j)
END DO
END DO
IF (jj < ddapc_restraint_control%natoms .AND. jj /= 0) &
IF (jj < ddapc_restraint_control%natoms .AND. jj /= 0) THEN
CPABORT("Need no or the same number of coeff as there are atoms.")
END IF
END DO
k = 0
DO i = 1, SIZE(qs_control%ddapc_restraint_control)
@ -3365,12 +3409,13 @@ CONTAINS
proj_mo_section => section_vals_get_subs_vals(rtp_section, "PRINT%PROJECTION_MO")
CALL section_vals_get(proj_mo_section, explicit=is_present)
IF (is_present) THEN
IF (dft_control%rtp_control%linear_scaling) &
IF (dft_control%rtp_control%linear_scaling) THEN
CALL cp_abort(__LOCATION__, &
"You have defined a time dependent projection of mos, but "// &
"only the density matrix is propagated (DENSITY_PROPAGATION "// &
".TRUE.). Please either use MO-based real time DFT or do not "// &
"define any PRINT%PROJECTION_MO section")
END IF
dft_control%rtp_control%is_proj_mo = .TRUE.
ELSE
dft_control%rtp_control%is_proj_mo = .FALSE.
@ -3455,8 +3500,9 @@ CONTAINS
DO i = 1, n_elems
DO j = 1, 2
IF (dft_control%rtp_control%print_pol_elements(i, j) > 3 .OR. &
dft_control%rtp_control%print_pol_elements(i, j) < 1) &
dft_control%rtp_control%print_pol_elements(i, j) < 1) THEN
CPABORT("Polarisation tensor element not 1,2 or 3 in at least one index")
END IF
END DO
END DO
END IF
@ -3579,8 +3625,9 @@ CONTAINS
END DO
dft_control%probe(i)%natoms = jj
IF (dft_control%probe(i)%natoms < 1) &
IF (dft_control%probe(i)%natoms < 1) THEN
CPABORT("Need at least 1 atom to use hair probes formalism")
END IF
ALLOCATE (dft_control%probe(i)%atom_ids(dft_control%probe(i)%natoms))
jj = 0

3
src/cp_dbcsr_cholesky.F
View file

@ -214,8 +214,9 @@ CONTAINS
CALL copy_dbcsr_to_fm(matrixb, fm_matrixb)
!CALL copy_dbcsr_to_fm(matrixout, fm_matrixout)
IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
IF (op /= "SOLVE" .AND. op /= "MULTIPLY") THEN
CPABORT("wrong argument op")
END IF
IF (PRESENT(pos)) THEN
SELECT CASE (pos)

46
src/cp_dbcsr_operations.F
View file

@ -611,8 +611,9 @@ CONTAINS
row_blk_size, col_blk_size_right_out)
CALL copy_fm_to_dbcsr(fm_in, in)
IF (ncol /= k_out .OR. my_beta /= 0.0_dp) &
IF (ncol /= k_out .OR. my_beta /= 0.0_dp) THEN
CALL copy_fm_to_dbcsr(fm_out, out)
END IF
CALL timeset(routineN//'_core', timing_handle_mult)
CALL dbcsr_multiply("N", "N", my_alpha, matrix, in, my_beta, out, &
@ -648,8 +649,9 @@ CONTAINS
n1 = SIZE(sizes1)
n2 = SIZE(sizes2)
IF (n1 /= n2) &
IF (n1 /= n2) THEN
CPABORT("distributions must be equal!")
END IF
sizes1(1:n1) = sizes2(1:n1)
used = SUM(sizes1(1:n1))
! If sizes1 does not cover everything, then we increase the
@ -723,8 +725,9 @@ CONTAINS
NULLIFY (col_dist_left)
IF (ncol > 0) THEN
IF (.NOT. dbcsr_valid_index(sparse_matrix)) &
IF (.NOT. dbcsr_valid_index(sparse_matrix)) THEN
CPABORT("sparse_matrix must pre-exist")
END IF
!
! Setup matrix_v
CALL cp_fm_get_info(matrix_v, ncol_global=k)
@ -830,25 +833,6 @@ CONTAINS
WRITE (*, *) 'PRESENT (matrix_g)', PRESENT(matrix_g)
WRITE (*, *) 'matrix_type=', dbcsr_get_matrix_type(sparse_matrix)
WRITE (*, *) 'norm(sm+alpha*v*g^t - fm+alpha*v*g^t)/n=', norm/REAL(nao, dp)
IF (norm/REAL(nao, dp) > 1e-12_dp) THEN
!WRITE(*,*) 'fm_matrix'
!DO j=1,SIZE(fm_matrix%local_data,2)
! DO i=1,SIZE(fm_matrix%local_data,1)
! WRITE(*,'(A,I3,A,I3,A,E26.16,A)') 'a(',i,',',j,')=',fm_matrix%local_data(i,j),';'
! ENDDO
!ENDDO
!WRITE(*,*) 'mat_v'
!CALL dbcsr_print(mat_v)
!WRITE(*,*) 'mat_g'
!CALL dbcsr_print(mat_g)
!WRITE(*,*) 'sparse_matrix'
!CALL dbcsr_print(sparse_matrix)
!WRITE(*,*) 'sparse_matrix2 (-sm + sparse(fm))'
!CALL dbcsr_print(sparse_matrix2)
!WRITE(*,*) 'sparse_matrix3 (copy of sm input)'
!CALL dbcsr_print(sparse_matrix3)
!stop
END IF
CALL dbcsr_release(sparse_matrix2)
CALL dbcsr_release(sparse_matrix3)
CALL cp_fm_release(fm_matrix)
@ -1025,11 +1009,13 @@ CONTAINS
estimated_blocks = max_blocks_per_bin*nbins
ALLOCATE (blk_dist(estimated_blocks), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_dist")
END IF
ALLOCATE (blk_sizes(estimated_blocks), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_sizes")
END IF
element_stack = 0
nblks = 0
DO blk_layer = 1, max_blocks_per_bin
@ -1051,23 +1037,27 @@ CONTAINS
block_size => blk_sizes
ELSE
ALLOCATE (block_distribution(nblks), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_dist")
END IF
block_distribution(:) = blk_dist(1:nblks)
DEALLOCATE (blk_dist)
ALLOCATE (block_size(nblks), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_sizes")
END IF
block_size(:) = blk_sizes(1:nblks)
DEALLOCATE (blk_sizes)
END IF
ELSE
ALLOCATE (block_distribution(0), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_dist")
END IF
ALLOCATE (block_size(0), stat=stat)
IF (stat /= 0) &
IF (stat /= 0) THEN
CPABORT("blk_sizes")
END IF
END IF
1579 FORMAT(I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5, 1X, I5)
IF (debug_mod) THEN

2
src/cp_ddapc_forces.F
View file

@ -217,7 +217,7 @@ CONTAINS
END IF
IF (iparticle1 /= iparticle2) THEN
ra = rvec
r = SQRT(DOT_PRODUCT(ra, ra))
r = NORM2(ra)
t2 = -1.0_dp/(r*r)*factor
drvec = ra/r*q1t*q2t
d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + t2*drvec

4
src/cp_ddapc_methods.F
View file

@ -204,7 +204,7 @@ CONTAINS
iparticle2, istart_g, s_dim
REAL(KIND=dp) :: g2, gcut2, tmp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_Am, my_Amw
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gfunc_sq(:, :, :)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gfunc_sq
!NB precalculate as many things outside of the innermost loop as possible, in particular w(ig)*gfunc(ig,igauus1)*gfunc(ig,igauss2)
@ -766,7 +766,7 @@ CONTAINS
!
IF (iparticle1 /= iparticle2) THEN
ra = rvec
r = SQRT(DOT_PRODUCT(ra, ra))
r = NORM2(ra)
my_val = factor/r
END IF
EwM(idim) = my_val - factor*g_ewald

9
src/cp_ddapc_util.F
View file

@ -355,7 +355,8 @@ CONTAINS
bv(:) = 0.0_dp
cv(:) = 1.0_dp/Vol
CALL build_b_vector(bv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
particle_set, radii, rho_tot_g, gcut); bv(:) = bv(:)/Vol
particle_set, radii, rho_tot_g, gcut)
bv(:) = bv(:)/Vol
CALL rho_tot_g%pw_grid%para%group%sum(bv)
c1 = DOT_PRODUCT(cv, MATMUL(cp_ddapc_env%AmI, bv)) - ch_dens
c1 = c1/cp_ddapc_env%c0
@ -718,11 +719,13 @@ CONTAINS
bv2(:) = 0.0_dp
particle_set(iparticle)%r(i) = rvec(i) + dx
CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
particle_set, radii, rho_tot_g, gcut); bv1(:) = bv1(:)/Vol
particle_set, radii, rho_tot_g, gcut)
bv1(:) = bv1(:)/Vol
CALL rho_tot_g%pw_grid%para%group%sum(bv1)
particle_set(iparticle)%r(i) = rvec(i) - dx
CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
particle_set, radii, rho_tot_g, gcut); bv2(:) = bv2(:)/Vol
particle_set, radii, rho_tot_g, gcut)
bv2(:) = bv2(:)/Vol
CALL rho_tot_g%pw_grid%para%group%sum(bv2)
ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx)
DO kk = 1, SIZE(ddbv)

6
src/cp_eri_mme_interface.F
View file

@ -563,8 +563,9 @@ CONTAINS
WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Cutoff error:", param%par%err_c
WRITE (unit_nr, '(T2, A, T72, ES9.2)') "ERI_MME| Total error (minimax + cutoff):", param%par%err_mm + param%par%err_c
END IF
IF (param%par%print_calib) &
IF (param%par%print_calib) THEN
WRITE (unit_nr, '(T2, A, T68, F13.10)') "ERI_MME| Minimax scaling constant in AM-GM estimate:", param%par%C_mm
END IF
END IF
END IF
@ -786,8 +787,9 @@ CONTAINS
CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL", i_val=potential)
CALL section_vals_val_get(eri_mme_test_section, "POTENTIAL_PARAM", r_val=pot_par)
IF (nzet <= 0) &
IF (nzet <= 0) THEN
CPABORT("Number of exponents NZET must be greater than 0.")
END IF
CALL init_orbital_pointers(l_max)

5
src/cp_external_control.F
View file

@ -63,8 +63,9 @@ CONTAINS
external_comm = comm
external_master_id = in_external_master_id
IF (PRESENT(in_scf_energy_message_tag)) &
IF (PRESENT(in_scf_energy_message_tag)) THEN
scf_energy_message_tag = in_scf_energy_message_tag
END IF
IF (PRESENT(in_exit_tag)) THEN
! the exit tag should be different from the mpi_probe tag default
CPASSERT(in_exit_tag /= -1)
@ -203,7 +204,7 @@ CONTAINS
IF (PRESENT(target_time)) THEN
my_target_time = target_time
my_start_time = start_time
ELSEIF (PRESENT(globenv)) THEN
ELSE IF (PRESENT(globenv)) THEN
my_target_time = globenv%cp2k_target_time
my_start_time = globenv%cp2k_start_time
ELSE

12
src/cp_subsys_methods.F
View file

@ -128,16 +128,19 @@ CONTAINS
subsys%para_env => para_env
my_use_motion_section = .FALSE.
IF (PRESENT(use_motion_section)) &
IF (PRESENT(use_motion_section)) THEN
my_use_motion_section = use_motion_section
END IF
my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
IF (PRESENT(force_env_section)) &
IF (PRESENT(force_env_section)) THEN
my_force_env_section => force_env_section
END IF
my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
IF (PRESENT(subsys_section)) &
IF (PRESENT(subsys_section)) THEN
my_subsys_section => subsys_section
END IF
CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
IF (SIZE(seed_vals) == 1) THEN
@ -288,8 +291,9 @@ CONTAINS
CPASSERT(.NOT. ASSOCIATED(small_subsys))
CPASSERT(ASSOCIATED(big_subsys))
IF (big_subsys%para_env /= small_para_env) &
IF (big_subsys%para_env /= small_para_env) THEN
CPABORT("big_subsys%para_env==small_para_env")
END IF
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------

2
src/cryssym.F
View file

@ -1755,7 +1755,7 @@ CONTAINS
csym%kplink(1, j) = i
wkp(i) = wkp(i) + 1.0_dp
kpop(j) = kr
ELSEIF (csym%kplink(1, j) /= i) THEN
ELSE IF (csym%kplink(1, j) /= i) THEN
! Approximate K290 operation sets need not be closed for structures whose
! coordinates lie close to several symmetry tolerances. Keep the existing
! disjoint orbit instead of aborting or double-counting this mesh point.

3
src/csvr_system_types.F
View file

@ -148,8 +148,9 @@ CONTAINS
SUBROUTINE csvr_thermo_dealloc(nvt)
TYPE(csvr_thermo_type), DIMENSION(:), POINTER :: nvt
IF (ASSOCIATED(nvt)) &
IF (ASSOCIATED(nvt)) THEN
DEALLOCATE (nvt)
END IF
END SUBROUTINE csvr_thermo_dealloc
END MODULE csvr_system_types

10
src/ct_methods.F
View file

@ -39,6 +39,7 @@ MODULE ct_methods
USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
USE kinds, ONLY: dp
USE machine, ONLY: m_walltime
USE mathconstants, ONLY: pi
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -1378,14 +1379,12 @@ CONTAINS
INTEGER, INTENT(OUT) :: nmins
INTEGER :: i, nroots
REAL(KIND=dp) :: DD, der, p, phi, pi, q, temp1, temp2, u, &
v, y1, y2, y2i, y2r, y3
REAL(KIND=dp) :: DD, der, p, phi, q, temp1, temp2, u, v, &
y1, y2, y2i, y2r, y3
REAL(KIND=dp), DIMENSION(3) :: x
! CALL timeset(routineN,handle)
pi = ACOS(-1.0_dp)
! Step 0: Check coefficients and find the true order of the eq
IF (a == 0.0_dp) THEN
IF (b == 0.0_dp) THEN
@ -1520,8 +1519,9 @@ CONTAINS
! create a matrix for eigenvectors
CALL dbcsr_work_create(c, work_mutable=.TRUE.)
IF (do_eigenvalues) &
IF (do_eigenvalues) THEN
CALL dbcsr_work_create(e, work_mutable=.TRUE.)
END IF
CALL dbcsr_iterator_readonly_start(iter, matrix)

68
src/ct_types.F
View file

@ -66,31 +66,6 @@ MODULE ct_types
REAL(KIND=dp) :: energy_correction = 0.0_dp
!SPIN!!! ! metric matrices for covariant to contravariant transformations
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: p_index_up=>NULL()
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: p_index_down=>NULL()
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: q_index_up=>NULL()
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: q_index_down=>NULL()
!SPIN!!!
!SPIN!!! ! kohn-sham, covariant-covariant representation
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_ks=>NULL()
!SPIN!!! ! density, contravariant-contravariant representation
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_p=>NULL()
!SPIN!!! ! occ orbitals, contravariant-covariant representation
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_t=>NULL()
!SPIN!!! ! virt orbitals, contravariant-covariant representation
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_v=>NULL()
!SPIN!!!
!SPIN!!! ! to avoid building Occ-by-N and Virt-vy-N matrices inside
!SPIN!!! ! the ct routines get them from the external code
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_qp_template=>NULL()
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER :: matrix_pq_template=>NULL()
!SPIN!!!
!SPIN!!! ! single excitation amplitudes
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), ALLOCATABLE :: matrix_x
!SPIN!!! ! residuals
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), ALLOCATABLE :: matrix_res
! metric matrices for covariant to contravariant transformations
TYPE(dbcsr_type), POINTER :: p_index_up => NULL()
TYPE(dbcsr_type), POINTER :: p_index_down => NULL()
@ -154,7 +129,6 @@ CONTAINS
env%order_lanczos = 3
env%eps_lancsoz = 1.0E-4_dp
env%max_iter_lanczos = 40
!env%nspins = -1
env%converged = .FALSE.
env%conjugator = cg_polak_ribiere
@ -231,22 +205,6 @@ CONTAINS
qq_preconditioner_full, &
pp_preconditioner_full
!INTEGER , OPTIONAL :: nspins
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: p_index_up
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: p_index_down
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: q_index_up
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: q_index_down
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_ks
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_p
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_t
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_v
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_qp_template
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_pq_template
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), POINTER, OPTIONAL :: matrix_x
!SPIN!!!
!SPIN!!! TYPE(dbcsr_type), DIMENSION(:), OPTIONAL :: copy_matrix_x
!INTEGER :: ispin
IF (PRESENT(use_occ_orbs)) use_occ_orbs = env%use_occ_orbs
IF (PRESENT(use_virt_orbs)) use_virt_orbs = env%use_virt_orbs
IF (PRESENT(occ_orbs_orthogonal)) occ_orbs_orthogonal = &
@ -267,7 +225,6 @@ CONTAINS
IF (PRESENT(eps_convergence)) eps_convergence = env%eps_convergence
IF (PRESENT(eps_filter)) eps_filter = env%eps_filter
IF (PRESENT(max_iter)) max_iter = env%max_iter
!IF (PRESENT(nspins)) nspins = env%nspins
IF (PRESENT(matrix_ks)) matrix_ks => env%matrix_ks
IF (PRESENT(matrix_p)) matrix_p => env%matrix_p
IF (PRESENT(matrix_t)) matrix_t => env%matrix_t
@ -281,12 +238,8 @@ CONTAINS
IF (PRESENT(p_index_down)) p_index_down => env%p_index_down
IF (PRESENT(q_index_down)) q_index_down => env%q_index_down
IF (PRESENT(copy_matrix_x)) THEN
!DO ispin=1,env%nspins
!CALL dbcsr_copy(copy_matrix_x(ispin),env%matrix_x(ispin))
CALL dbcsr_copy(copy_matrix_x, env%matrix_x)
!ENDDO
END IF
!IF (PRESENT(matrix_x)) matrix_x => env%matrix_x
IF (PRESENT(energy_correction)) energy_correction = env%energy_correction
IF (PRESENT(converged)) converged = env%converged
@ -352,20 +305,6 @@ CONTAINS
LOGICAL, OPTIONAL :: qq_preconditioner_full, &
pp_preconditioner_full
!INTEGER , OPTIONAL :: nspins
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: p_index_up
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: p_index_down
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: q_index_up
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: q_index_down
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_ks
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_p
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_t
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_v
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_qp_template
!SPIN!!! TYPE(dbcsr_type), TARGET, DIMENSION(:), OPTIONAL :: matrix_pq_template
! set para_env and blacs_env which are needed to operate with full matrices
! it would be nice to have everything with cp_dbcsr matrices, well maybe later
env%para_env => para_env
env%blacs_env => blacs_env
@ -389,7 +328,6 @@ CONTAINS
IF (PRESENT(eps_convergence)) env%eps_convergence = eps_convergence
IF (PRESENT(eps_filter)) env%eps_filter = eps_filter
IF (PRESENT(max_iter)) env%max_iter = max_iter
!IF (PRESENT(nspins)) env%nspins = nspins
IF (PRESENT(conjugator)) env%conjugator = conjugator
IF (PRESENT(matrix_ks)) env%matrix_ks => matrix_ks
IF (PRESENT(matrix_p)) env%matrix_p => matrix_p
@ -415,18 +353,12 @@ CONTAINS
TYPE(ct_step_env_type) :: env
!INTEGER :: ispin
NULLIFY (env%para_env)
NULLIFY (env%blacs_env)
!DO ispin=1,env%nspins
CALL dbcsr_release(env%matrix_x)
CALL dbcsr_release(env%matrix_res)
!CALL dbcsr_release(env%matrix_x(ispin))
!CALL dbcsr_release(env%matrix_res(ispin))
!ENDDO
!DEALLOCATE(env%matrix_x,env%matrix_res)
NULLIFY (env%p_index_up)
NULLIFY (env%p_index_down)

3
src/dbm/dbm_api.F
View file

@ -182,8 +182,9 @@ CONTAINS
num_blocks_diff = ABS(num_blocks - num_blocks_dbcsr)
IF (num_blocks_diff /= 0) THEN
WRITE (*, *) "num_blocks mismatch dbcsr:", num_blocks_dbcsr, "new:", num_blocks
IF (DBM_VALIDATE_NBLOCKS_MATCH) &
IF (DBM_VALIDATE_NBLOCKS_MATCH) THEN
CPABORT("num_blocks mismatch")
END IF
END IF
IF (DBM_VALIDATE_NBLOCKS_MATCH) THEN

2
src/dbm/dbm_tests.F
View file

@ -439,7 +439,7 @@ CONTAINS
INTEGER(KIND=int_8) :: map
map = ((irow - 1 + icol*INT(nrow, int_8))*(1 + MODULO(ival, 2**16)))*2 + 1 + 0*ncol ! ncol used
iseed(4) = INT(MODULO(map, 2_int_8**12)); map = map/2_int_8**12; ! keep odd
iseed(4) = INT(MODULO(map, 2_int_8**12)); map = map/2_int_8**12 ! keep odd
iseed(3) = INT(MODULO(IEOR(map, 3541_int_8), 2_int_8**12)); map = map/2_int_8**12
iseed(2) = INT(MODULO(IEOR(map, 1153_int_8), 2_int_8**12)); map = map/2_int_8**12
iseed(1) = INT(MODULO(IEOR(map, 2029_int_8), 2_int_8**12)); map = map/2_int_8**12

2
src/dbt/dbt_allocate_wrap.F
View file

@ -56,7 +56,7 @@ CONTAINS
ELSE
shape_prv = shape_spec
END IF
ELSEIF (PRESENT(source)) THEN
ELSE IF (PRESENT(source)) THEN
IF (PRESENT(order)) THEN
shape_prv(order) = SHAPE(source)
ELSE

22
src/dbt/dbt_array_list_methods.F
View file

@ -51,7 +51,7 @@ MODULE dbt_array_list_methods
TYPE array_list
INTEGER, DIMENSION(:), ALLOCATABLE :: col_data
INTEGER, DIMENSION(:), ALLOCATABLE :: ptr
END TYPE
END TYPE array_list
INTERFACE get_ith_array
MODULE PROCEDURE allocate_and_get_ith_array
@ -128,7 +128,7 @@ CONTAINS
END IF
#:endfor
END SUBROUTINE
END SUBROUTINE create_array_list
! **************************************************************************************************
!> \brief extract a subset of arrays
@ -151,7 +151,7 @@ CONTAINS
CALL create_array_list(array_sublist, ndata, ${varlist("data", nmax=dim)}$)
END IF
#:endfor
END FUNCTION
END FUNCTION array_sublist
! **************************************************************************************************
!> \brief destroy array list.
@ -161,7 +161,7 @@ CONTAINS
TYPE(array_list), INTENT(INOUT) :: list
DEALLOCATE (list%ptr, list%col_data)
END SUBROUTINE
END SUBROUTINE destroy_array_list
! **************************************************************************************************
!> \brief Get all arrays contained in list
@ -185,7 +185,7 @@ CONTAINS
o(1:ndata) = i_selected(:)
ELSE
ndata = number_of_arrays(list)
o(1:ndata) = (/(i, i=1, ndata)/)
o(1:ndata) = [(i, i=1, ndata)]
END IF
ASSOCIATE (ptr => list%ptr, col_data => list%col_data)
@ -215,7 +215,7 @@ CONTAINS
END ASSOCIATE
END SUBROUTINE
END SUBROUTINE get_ith_array
! **************************************************************************************************
!> \brief get ith array
@ -231,7 +231,7 @@ CONTAINS
ALLOCATE (array, source=col_data(ptr(i):ptr(i + 1) - 1))
END ASSOCIATE
END SUBROUTINE
END SUBROUTINE allocate_and_get_ith_array
! **************************************************************************************************
!> \brief sizes of arrays stored in list
@ -288,7 +288,7 @@ CONTAINS
partial_sum = partial_sum + list_in%col_data(i_ptr)
END DO
END DO
END SUBROUTINE
END SUBROUTINE array_offsets
! **************************************************************************************************
!> \brief reorder array list.
@ -309,7 +309,7 @@ CONTAINS
END IF
#:endfor
END SUBROUTINE
END SUBROUTINE reorder_arrays
! **************************************************************************************************
!> \brief check whether two array lists are equal
@ -320,7 +320,7 @@ CONTAINS
LOGICAL :: check_equal
check_equal = array_eq_i(list1%col_data, list2%col_data) .AND. array_eq_i(list1%ptr, list2%ptr)
END FUNCTION
END FUNCTION check_equal
! **************************************************************************************************
!> \brief check whether two arrays are equal
@ -337,6 +337,6 @@ CONTAINS
array_eq_i = .FALSE.
IF (SIZE(arr1) == SIZE(arr2)) array_eq_i = ALL(arr1 == arr2)
#endif
END FUNCTION
END FUNCTION array_eq_i
END MODULE dbt_array_list_methods

16
src/dbt/dbt_methods.F
View file

@ -231,7 +231,7 @@ CONTAINS
IF (.NOT. PRESENT(order)) THEN
IF (array_eq_i(map1_in_1, map2_in_1) .AND. array_eq_i(map1_in_2, map2_in_2)) THEN
dist_compatible_tas = check_equal(in_tmp_3%nd_dist, out_tmp_1%nd_dist)
ELSEIF (array_eq_i([map1_in_1, map1_in_2], [map2_in_1, map2_in_2])) THEN
ELSE IF (array_eq_i([map1_in_1, map1_in_2], [map2_in_1, map2_in_2])) THEN
dist_compatible_tensor = check_equal(in_tmp_3%nd_dist, out_tmp_1%nd_dist)
END IF
END IF
@ -239,7 +239,7 @@ CONTAINS
IF (dist_compatible_tas) THEN
CALL dbt_tas_copy(out_tmp_1%matrix_rep, in_tmp_3%matrix_rep, summation)
IF (move_prv) CALL dbt_clear(in_tmp_3)
ELSEIF (dist_compatible_tensor) THEN
ELSE IF (dist_compatible_tensor) THEN
CALL dbt_copy_nocomm(in_tmp_3, out_tmp_1, summation)
IF (move_prv) CALL dbt_clear(in_tmp_3)
ELSE
@ -1274,7 +1274,7 @@ CONTAINS
ALLOCATE (tensor2_out)
CALL dbt_remap(tensor2, ind2_linked, ind2_free, tensor2_out, comm_2d=dist_in%pgrid%mp_comm_2d, &
dist1=dist_list, mp_dims_1=mp_dims, nodata=nodata2, move_data=move_data_2)
ELSEIF (compat1 == 2) THEN ! linked index is second 2d dimension
ELSE IF (compat1 == 2) THEN ! linked index is second 2d dimension
! get distribution of linked index, tensor 2 must adopt this distribution
! get grid dimensions of linked index
ALLOCATE (mp_dims(ndims_mapping_column(dist_in%pgrid%nd_index_grid)))
@ -1315,7 +1315,7 @@ CONTAINS
ALLOCATE (tensor1_out)
CALL dbt_remap(tensor1, ind1_linked, ind1_free, tensor1_out, comm_2d=dist_in%pgrid%mp_comm_2d, &
dist1=dist_list, mp_dims_1=mp_dims, nodata=nodata1, move_data=move_data_1)
ELSEIF (compat2 == 2) THEN
ELSE IF (compat2 == 2) THEN
ALLOCATE (mp_dims(ndims_mapping_column(dist_in%pgrid%nd_index_grid)))
CALL dbt_get_mapping_info(dist_in%pgrid%nd_index_grid, dims2_2d=mp_dims)
ALLOCATE (tensor1_out)
@ -1419,7 +1419,7 @@ CONTAINS
WRITE (unit_nr_prv, '(T2,A,1X,A,A,1X)', advance='no') "compatibility of", TRIM(tensor_in%name), ":"
IF (compat1 == 1 .AND. compat2 == 2) THEN
WRITE (unit_nr_prv, '(A)') "Normal"
ELSEIF (compat1 == 2 .AND. compat2 == 1) THEN
ELSE IF (compat1 == 2 .AND. compat2 == 1) THEN
WRITE (unit_nr_prv, '(A)') "Transposed"
ELSE
WRITE (unit_nr_prv, '(A)') "Not compatible"
@ -1442,7 +1442,7 @@ CONTAINS
IF (compat1 == 1 .AND. compat2 == 2) THEN
trans = .FALSE.
ELSEIF (compat1 == 2 .AND. compat2 == 1) THEN
ELSE IF (compat1 == 2 .AND. compat2 == 1) THEN
trans = .TRUE.
ELSE
CPABORT("this should not happen")
@ -1453,7 +1453,7 @@ CONTAINS
WRITE (unit_nr_prv, '(T2,A,1X,A,A,1X)', advance='no') "compatibility of", TRIM(tensor_out%name), ":"
IF (compat1 == 1 .AND. compat2 == 2) THEN
WRITE (unit_nr_prv, '(A)') "Normal"
ELSEIF (compat1 == 2 .AND. compat2 == 1) THEN
ELSE IF (compat1 == 2 .AND. compat2 == 1) THEN
WRITE (unit_nr_prv, '(A)') "Transposed"
ELSE
WRITE (unit_nr_prv, '(A)') "Not compatible"
@ -1526,7 +1526,7 @@ CONTAINS
compat_map = 0
IF (array_eq_i(map1, compat_ind)) THEN
compat_map = 1
ELSEIF (array_eq_i(map2, compat_ind)) THEN
ELSE IF (array_eq_i(map2, compat_ind)) THEN
compat_map = 2
END IF

17
src/dbt/dbt_split.F
View file

@ -272,7 +272,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE dbt_split_blocks_generic
! **************************************************************************************************
!> \brief Split tensor blocks into smaller blocks of maximum size PRODUCT(block_sizes).
@ -331,7 +331,7 @@ CONTAINS
nodata=nodata)
#:endfor
END SUBROUTINE
END SUBROUTINE dbt_split_blocks
! **************************************************************************************************
!> \brief Copy tensor with split blocks to tensor with original block sizes.
@ -480,7 +480,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE dbt_split_copyback
! **************************************************************************************************
!> \brief split two tensors with same total sizes but different block sizes such that they have
@ -503,7 +503,8 @@ CONTAINS
INTEGER, DIMENSION(ndims_tensor(tensor1)), &
INTENT(IN), OPTIONAL :: order
LOGICAL, INTENT(IN), OPTIONAL :: nodata1, nodata2, move_data
INTEGER, DIMENSION(:), ALLOCATABLE :: ${varlist("blk_size_split_1")}$, ${varlist("blk_size_split_2")}$, &
INTEGER, DIMENSION(:), ALLOCATABLE :: ${varlist("blk_size_split_1")}$, &
${varlist("blk_size_split_2")}$, &
blk_size_d_1, blk_size_d_2, blk_size_d_split
INTEGER :: size_sum_1, size_sum_2, size_sum, bind_1, bind_2, isplit, bs, idim, i
LOGICAL :: move_prv, nodata1_prv, nodata2_prv
@ -551,7 +552,7 @@ CONTAINS
size_sum = size_sum + bs
isplit = isplit + 1
blk_size_d_split(isplit) = bs
ELSEIF (blk_size_d_1(bind_1 + 1) > blk_size_d_2(bind_2 + 1)) THEN
ELSE IF (blk_size_d_1(bind_1 + 1) > blk_size_d_2(bind_2 + 1)) THEN
bind_2 = bind_2 + 1
bs = blk_size_d_2(bind_2)
blk_size_d_1(bind_1 + 1) = blk_size_d_1(bind_1 + 1) - bs
@ -608,7 +609,7 @@ CONTAINS
END IF
#:endfor
END SUBROUTINE
END SUBROUTINE dbt_make_compatible_blocks
! **************************************************************************************************
!> \author Patrick Seewald
@ -698,6 +699,6 @@ CONTAINS
CALL dbt_copy_contraction_storage(tensor_in, tensor_out)
CALL timestop(handle)
END SUBROUTINE
END SUBROUTINE dbt_crop
END MODULE
END MODULE dbt_split

10
src/dbt/dbt_types.F
View file

@ -211,7 +211,7 @@ CONTAINS
CALL dbt_get_mapping_info(map_grid, &
dims_2d=grid_dims, &
dims1_2d=new_dbt_tas_dist_t%dims_grid)
ELSEIF (which_dim == 2) THEN
ELSE IF (which_dim == 2) THEN
ALLOCATE (new_dbt_tas_dist_t%dims(ndims_mapping_column(map_blks)))
ALLOCATE (index_map(ndims_mapping_column(map_blks)))
CALL dbt_get_mapping_info(map_blks, &
@ -335,7 +335,7 @@ CONTAINS
dims_2d_i8=matrix_dims, &
map1_2d=index_map, &
dims1_2d=new_dbt_tas_blk_size_t%dims)
ELSEIF (which_dim == 2) THEN
ELSE IF (which_dim == 2) THEN
ALLOCATE (index_map(ndims_mapping_column(map_blks)))
ALLOCATE (new_dbt_tas_blk_size_t%dims(ndims_mapping_column(map_blks)))
CALL dbt_get_mapping_info(map_blks, &
@ -455,7 +455,7 @@ CONTAINS
IF (idim /= SIZE(tensor_dims_sorted)) THEN
dims(idim + 1:) = 0
CALL mp_dims_create(pdims_rem, dims(idim + 1:))
ELSEIF (lb_ratio_prv < 0.5_dp) THEN
ELSE IF (lb_ratio_prv < 0.5_dp) THEN
! resort to a less strict load imbalance factor
dims(:) = dims_store
CALL dbt_mp_dims_create(nodes, dims, tensor_dims, 0.5_dp)
@ -836,7 +836,7 @@ CONTAINS
abort = .FALSE.
IF (.NOT. ASSOCIATED(dist%refcount)) THEN
abort = .TRUE.
ELSEIF (dist%refcount < 1) THEN
ELSE IF (dist%refcount < 1) THEN
abort = .TRUE.
END IF
@ -1264,7 +1264,7 @@ CONTAINS
abort = .FALSE.
IF (.NOT. ASSOCIATED(tensor%refcount)) THEN
abort = .TRUE.
ELSEIF (tensor%refcount < 1) THEN
ELSE IF (tensor%refcount < 1) THEN
abort = .TRUE.
END IF

16
src/dbt/tas/dbt_tas_mm.F
View file

@ -276,7 +276,7 @@ CONTAINS
WRITE (unit_nr_prv, "(T4,A,T68,I13)") "Est. optimal split factor:", nsplit
END IF
ELSEIF (batched_repl > 0) THEN
ELSE IF (batched_repl > 0) THEN
nsplit = nsplit_batched
nsplit_opt = nsplit
max_mm_dim = max_mm_dim_batched
@ -342,7 +342,7 @@ CONTAINS
IF (matrix_c%do_batched == 1) THEN
matrix_c%mm_storage%batched_beta = beta
ELSEIF (matrix_c%do_batched > 1) THEN
ELSE IF (matrix_c%do_batched > 1) THEN
matrix_c%mm_storage%batched_beta = matrix_c%mm_storage%batched_beta*beta
END IF
@ -355,7 +355,7 @@ CONTAINS
IF (.NOT. nodata_3) CALL dbm_zero(matrix_c_rs%matrix)
IF (matrix_c%do_batched >= 1) matrix_c%mm_storage%store_batched => matrix_c_rs
ELSEIF (matrix_c%do_batched == 3) THEN
ELSE IF (matrix_c%do_batched == 3) THEN
matrix_c_rs => matrix_c%mm_storage%store_batched
END IF
@ -446,7 +446,7 @@ CONTAINS
matrix_b%mm_storage%store_batched_repl => matrix_b_rep
CALL dbt_tas_set_batched_state(matrix_b, state=3)
END IF
ELSEIF (matrix_b%do_batched == 3) THEN
ELSE IF (matrix_b%do_batched == 3) THEN
matrix_b_rep => matrix_b%mm_storage%store_batched_repl
END IF
@ -532,14 +532,14 @@ CONTAINS
matrix_c%mm_storage%store_batched_repl => matrix_c_rep
CALL dbt_tas_set_batched_state(matrix_c, state=3)
END IF
ELSEIF (matrix_c%do_batched == 2) THEN
ELSE IF (matrix_c%do_batched == 2) THEN
ALLOCATE (matrix_c_rep)
CALL dbt_tas_replicate(matrix_c_rs%matrix, dbt_tas_info(matrix_a_rs), matrix_c_rep, nodata=nodata_3)
! just leave sparsity structure for retain sparsity but no values
IF (.NOT. nodata_3) CALL dbm_zero(matrix_c_rep%matrix)
matrix_c%mm_storage%store_batched_repl => matrix_c_rep
CALL dbt_tas_set_batched_state(matrix_c, state=3)
ELSEIF (matrix_c%do_batched == 3) THEN
ELSE IF (matrix_c%do_batched == 3) THEN
matrix_c_rep => matrix_c%mm_storage%store_batched_repl
END IF
@ -625,7 +625,7 @@ CONTAINS
matrix_a%mm_storage%store_batched_repl => matrix_a_rep
CALL dbt_tas_set_batched_state(matrix_a, state=3)
END IF
ELSEIF (matrix_a%do_batched == 3) THEN
ELSE IF (matrix_a%do_batched == 3) THEN
matrix_a_rep => matrix_a%mm_storage%store_batched_repl
END IF
@ -740,7 +740,7 @@ CONTAINS
CALL dbt_tas_destroy(matrix_c_rs)
DEALLOCATE (matrix_c_rs)
IF (PRESENT(filter_eps)) CALL dbt_tas_filter(matrix_c, filter_eps)
ELSEIF (matrix_c%do_batched > 0) THEN
ELSE IF (matrix_c%do_batched > 0) THEN
IF (matrix_c%mm_storage%batched_out) THEN
matrix_c%mm_storage%batched_trans = (transc_prv .NEQV. transc)
END IF

8
src/dbt/tas/dbt_tas_split.F
View file

@ -266,7 +266,7 @@ CONTAINS
nsplit_list_square(count_square) = split
count_accept = count_accept + 1
nsplit_list_accept(count_accept) = split
ELSEIF (accept_pgrid_dims(dims_sub, relative=.FALSE.)) THEN
ELSE IF (accept_pgrid_dims(dims_sub, relative=.FALSE.)) THEN
count_accept = count_accept + 1
nsplit_list_accept(count_accept) = split
END IF
@ -277,10 +277,10 @@ CONTAINS
IF (count_square > 0) THEN
minpos = MINLOC(ABS(nsplit_list_square(1:count_square) - nsplit), DIM=1)
get_opt_nsplit = nsplit_list_square(minpos)
ELSEIF (count_accept > 0) THEN
ELSE IF (count_accept > 0) THEN
minpos = MINLOC(ABS(nsplit_list_accept(1:count_accept) - nsplit), DIM=1)
get_opt_nsplit = nsplit_list_accept(minpos)
ELSEIF (count > 0) THEN
ELSE IF (count > 0) THEN
minpos = MINLOC(ABS(nsplit_list(1:count) - nsplit), DIM=1)
get_opt_nsplit = nsplit_list(minpos)
ELSE
@ -462,7 +462,7 @@ CONTAINS
IF (.NOT. ASSOCIATED(split_info%refcount)) THEN
abort = .TRUE.
ELSEIF (split_info%refcount < 1) THEN
ELSE IF (split_info%refcount < 1) THEN
abort = .TRUE.
END IF

8
src/dbt/tas/dbt_tas_util.F
View file

@ -53,7 +53,7 @@ CONTAINS
tmp = arr(1)
arr(1) = arr(2)
arr(2) = tmp
END SUBROUTINE
END SUBROUTINE swap_i8
! **************************************************************************************************
!> \brief ...
@ -68,7 +68,7 @@ CONTAINS
tmp = arr(1)
arr(1) = arr(2)
arr(2) = tmp
END SUBROUTINE
END SUBROUTINE swap_i
! **************************************************************************************************
!> \brief ...
@ -87,7 +87,7 @@ CONTAINS
array_eq_i = .FALSE.
IF (SIZE(arr1) == SIZE(arr2)) array_eq_i = ALL(arr1 == arr2)
#endif
END FUNCTION
END FUNCTION array_eq_i
! **************************************************************************************************
!> \brief ...
@ -106,6 +106,6 @@ CONTAINS
array_eq_i8 = .FALSE.
IF (SIZE(arr1) == SIZE(arr2)) array_eq_i8 = ALL(arr1 == arr2)
#endif
END FUNCTION
END FUNCTION array_eq_i8
END MODULE

3
src/dbx/cp_dbcsr_api.F
View file

@ -235,10 +235,11 @@ CONTAINS
TYPE(dbcsr_type), POINTER :: matrix
CALL dbcsr_release(matrix)
IF (dbcsr_valid_index(matrix)) &
IF (dbcsr_valid_index(matrix)) THEN
CALL cp_abort(__LOCATION__, &
'You should not "deallocate" a referenced matrix. '// &
'Avoid pointers to DBCSR matrices.')
END IF
DEALLOCATE (matrix)
END SUBROUTINE dbcsr_deallocate_matrix

10
src/debug_os_integrals.F
View file

@ -75,7 +75,7 @@ CONTAINS
ALLOCATE (swork(lds, lds, 1))
sab = 0._dp
rab(:) = B(:) - A(:)
dab = SQRT(DOT_PRODUCT(rab, rab))
dab = NORM2(rab)
xa_work(1) = xa
xb_work(1) = xb
rpgfa = 20._dp
@ -219,11 +219,11 @@ CONTAINS
!---------------------------------------
rab(:) = B(:) - A(:)
dab = SQRT(DOT_PRODUCT(rab, rab))
dab = NORM2(rab)
rac(:) = C(:) - A(:)
dac = SQRT(DOT_PRODUCT(rac, rac))
dac = NORM2(rac)
rbc(:) = C(:) - B(:)
dbc = SQRT(DOT_PRODUCT(rbc, rbc))
dbc = NORM2(rbc)
ALLOCATE (sabc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max)))
xa_work(1) = xa
xb_work(1) = xb
@ -422,7 +422,7 @@ CONTAINS
ALLOCATE (swork(lds, lds))
saabb = 0._dp
rab(:) = B(:) - A(:)
dab = SQRT(DOT_PRODUCT(rab, rab))
dab = NORM2(rab)
xa_work1(1) = xa1
xa_work2(1) = xa2
xb_work1(1) = xb1

8
src/dft_plus_u.F
View file

@ -189,7 +189,7 @@ CONTAINS
!> \date 02.07.2008
!> \par
!> \f{eqnarray*}{
!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}\\\
!> E^{\rm DFT+U} & = & E^{\rm DFT} + E^{\rm U}
!> & = & E^{\rm DFT} + \frac{1}{2}(U - J)\sum_\mu (q_\mu - q_\mu^2)\\[1ex]
!> V_{\mu\nu}^{\rm DFT+U} & = & V_{\mu\nu}^{\rm DFT} + V_{\mu\nu}^{\rm U}\\\
!> & = & \frac{\partial E^{\rm DFT}}
@ -766,10 +766,11 @@ CONTAINS
CALL para_env%sum(energy%dft_plus_u)
IF (energy%dft_plus_u < 0.0_dp) &
IF (energy%dft_plus_u < 0.0_dp) THEN
CALL cp_warn(__LOCATION__, &
"DFT+U energy contribution is negative possibly due "// &
"to unphysical Lowdin charges!")
END IF
! Release (local) full matrices
NULLIFY (fm_s_half)
@ -1991,10 +1992,11 @@ CONTAINS
CALL para_env%sum(energy%dft_plus_u)
IF (energy%dft_plus_u < 0.0_dp) &
IF (energy%dft_plus_u < 0.0_dp) THEN
CALL cp_warn(__LOCATION__, &
"DFT+U energy contribution is negative possibly due "// &
"to unphysical Mulliken charges!")
END IF
! Release local work storage

36
src/distribution_2d_types.F
View file

@ -132,8 +132,9 @@ CONTAINS
END IF
IF (PRESENT(n_col_distribution)) THEN
IF (ASSOCIATED(distribution_2d%col_distribution)) THEN
IF (n_col_distribution > distribution_2d%n_col_distribution) &
IF (n_col_distribution > distribution_2d%n_col_distribution) THEN
CPABORT("n_col_distribution<=distribution_2d%n_col_distribution")
END IF
! else alloc col_distribution?
END IF
distribution_2d%n_col_distribution = n_col_distribution
@ -145,15 +146,17 @@ CONTAINS
END IF
IF (PRESENT(n_row_distribution)) THEN
IF (ASSOCIATED(distribution_2d%row_distribution)) THEN
IF (n_row_distribution > distribution_2d%n_row_distribution) &
IF (n_row_distribution > distribution_2d%n_row_distribution) THEN
CPABORT("n_row_distribution<=distribution_2d%n_row_distribution")
END IF
! else alloc row_distribution?
END IF
distribution_2d%n_row_distribution = n_row_distribution
END IF
IF (PRESENT(local_rows_ptr)) &
IF (PRESENT(local_rows_ptr)) THEN
distribution_2d%local_rows => local_rows_ptr
END IF
IF (.NOT. ASSOCIATED(distribution_2d%local_rows)) THEN
CPASSERT(PRESENT(n_local_rows))
ALLOCATE (distribution_2d%local_rows(SIZE(n_local_rows)))
@ -164,11 +167,13 @@ CONTAINS
END IF
ALLOCATE (distribution_2d%n_local_rows(SIZE(distribution_2d%local_rows)))
IF (PRESENT(n_local_rows)) THEN
IF (SIZE(distribution_2d%n_local_rows) /= SIZE(n_local_rows)) &
IF (SIZE(distribution_2d%n_local_rows) /= SIZE(n_local_rows)) THEN
CPABORT("SIZE(distribution_2d%n_local_rows)==SIZE(n_local_rows)")
END IF
DO i = 1, SIZE(distribution_2d%n_local_rows)
IF (SIZE(distribution_2d%local_rows(i)%array) < n_local_rows(i)) &
IF (SIZE(distribution_2d%local_rows(i)%array) < n_local_rows(i)) THEN
CPABORT("SIZE(distribution_2d%local_rows(i)%array)>=n_local_rows(i)")
END IF
distribution_2d%n_local_rows(i) = n_local_rows(i)
END DO
ELSE
@ -178,8 +183,9 @@ CONTAINS
END DO
END IF
IF (PRESENT(local_cols_ptr)) &
IF (PRESENT(local_cols_ptr)) THEN
distribution_2d%local_cols => local_cols_ptr
END IF
IF (.NOT. ASSOCIATED(distribution_2d%local_cols)) THEN
CPASSERT(PRESENT(n_local_cols))
ALLOCATE (distribution_2d%local_cols(SIZE(n_local_cols)))
@ -190,11 +196,13 @@ CONTAINS
END IF
ALLOCATE (distribution_2d%n_local_cols(SIZE(distribution_2d%local_cols)))
IF (PRESENT(n_local_cols)) THEN
IF (SIZE(distribution_2d%n_local_cols) /= SIZE(n_local_cols)) &
IF (SIZE(distribution_2d%n_local_cols) /= SIZE(n_local_cols)) THEN
CPABORT("SIZE(distribution_2d%n_local_cols)==SIZE(n_local_cols)")
END IF
DO i = 1, SIZE(distribution_2d%n_local_cols)
IF (SIZE(distribution_2d%local_cols(i)%array) < n_local_cols(i)) &
IF (SIZE(distribution_2d%local_cols(i)%array) < n_local_cols(i)) THEN
CPABORT("SIZE(distribution_2d%local_cols(i)%array)>=n_local_cols(i)")
END IF
distribution_2d%n_local_cols(i) = n_local_cols(i)
END DO
ELSE
@ -315,8 +323,9 @@ CONTAINS
DO i = 1, SIZE(distribution_2d%row_distribution, 1)
WRITE (unit=unit_nr, fmt="(i6,',')", advance="no") distribution_2d%row_distribution(i, 1)
! keep lines finite, so that we can open outputs in vi
IF (MODULO(i, 8) == 0 .AND. i /= SIZE(distribution_2d%row_distribution, 1)) &
IF (MODULO(i, 8) == 0 .AND. i /= SIZE(distribution_2d%row_distribution, 1)) THEN
WRITE (unit=unit_nr, fmt='()')
END IF
END DO
WRITE (unit=unit_nr, fmt="('),')")
ELSE
@ -336,8 +345,9 @@ CONTAINS
DO i = 1, SIZE(distribution_2d%col_distribution, 1)
WRITE (unit=unit_nr, fmt="(i6,',')", advance="no") distribution_2d%col_distribution(i, 1)
! keep lines finite, so that we can open outputs in vi
IF (MODULO(i, 8) == 0 .AND. i /= SIZE(distribution_2d%col_distribution, 1)) &
IF (MODULO(i, 8) == 0 .AND. i /= SIZE(distribution_2d%col_distribution, 1)) THEN
WRITE (unit=unit_nr, fmt='()')
END IF
END DO
WRITE (unit=unit_nr, fmt="('),')")
ELSE
@ -355,8 +365,9 @@ CONTAINS
DO i = 1, SIZE(distribution_2d%n_local_rows)
WRITE (unit=unit_nr, fmt="(i6,',')", advance="no") distribution_2d%n_local_rows(i)
! keep lines finite, so that we can open outputs in vi
IF (MODULO(i, 10) == 0 .AND. i /= SIZE(distribution_2d%n_local_rows)) &
IF (MODULO(i, 10) == 0 .AND. i /= SIZE(distribution_2d%n_local_rows)) THEN
WRITE (unit=unit_nr, fmt='()')
END IF
END DO
WRITE (unit=unit_nr, fmt="('),')")
ELSE
@ -395,8 +406,9 @@ CONTAINS
DO i = 1, SIZE(distribution_2d%n_local_cols)
WRITE (unit=unit_nr, fmt="(i6,',')", advance="no") distribution_2d%n_local_cols(i)
! keep lines finite, so that we can open outputs in vi
IF (MODULO(i, 10) == 0 .AND. i /= SIZE(distribution_2d%n_local_cols)) &
IF (MODULO(i, 10) == 0 .AND. i /= SIZE(distribution_2d%n_local_cols)) THEN
WRITE (unit=unit_nr, fmt='()')
END IF
END DO
WRITE (unit=unit_nr, fmt="('),')")
ELSE

Some files were not shown because too many files have changed in this diff Show more