|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE s_square_methods 00013 00014 USE cp_control_types, ONLY: s2_restraint_type 00015 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 00016 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00017 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm 00018 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00019 cp_fm_struct_release,& 00020 cp_fm_struct_type 00021 USE cp_fm_types, ONLY: cp_fm_create,& 00022 cp_fm_get_info,& 00023 cp_fm_p_type,& 00024 cp_fm_release,& 00025 cp_fm_type 00026 USE cp_para_types, ONLY: cp_blacs_env_type,& 00027 cp_para_env_type 00028 USE f77_blas 00029 USE input_constants, ONLY: do_s2_constraint,& 00030 do_s2_restraint 00031 USE kinds, ONLY: dp 00032 USE message_passing, ONLY: mp_sum 00033 USE qs_mo_types, ONLY: get_mo_set,& 00034 mo_set_p_type 00035 USE termination, ONLY: stop_program 00036 USE timings, ONLY: timeset,& 00037 timestop 00038 #include "cp_common_uses.h" 00039 00040 IMPLICIT NONE 00041 00042 PRIVATE 00043 00044 ! *** Global parameters *** 00045 00046 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_square_methods' 00047 00048 PUBLIC :: compute_s_square,s2_restraint 00049 00050 CONTAINS 00051 00052 ! ***************************************************************************** 00064 SUBROUTINE compute_s_square(mos, matrix_s, s_square, s_square_ideal,& 00065 mo_derivs,strength,error) 00066 TYPE(mo_set_p_type), DIMENSION(:), 00067 POINTER :: mos 00068 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00069 POINTER :: matrix_s 00070 REAL(KIND=dp) :: s_square, s_square_ideal 00071 TYPE(cp_fm_p_type), DIMENSION(:), 00072 OPTIONAL, POINTER :: mo_derivs 00073 REAL(KIND=dp), OPTIONAL :: strength 00074 TYPE(cp_error_type), INTENT(INOUT) :: error 00075 00076 CHARACTER(len=*), PARAMETER :: routineN = 'compute_s_square', 00077 routineP = moduleN//':'//routineN 00078 00079 INTEGER :: handle, i, j, na, nalpha, nb, 00080 nbeta, ncol_local, nrow, 00081 nrow_local 00082 LOGICAL :: failure, uniform_occupation 00083 REAL(KIND=dp) :: tmp 00084 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 00085 TYPE(cp_blacs_env_type), POINTER :: context 00086 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00087 TYPE(cp_fm_type), POINTER :: c_alpha, c_beta, 00088 matrix_overlap, matrix_sc_a, 00089 matrix_sc_b 00090 TYPE(cp_para_env_type), POINTER :: para_env 00091 00092 CALL timeset(routineN,handle) 00093 00094 failure=.FALSE. 00095 00096 NULLIFY(fm_struct_tmp,matrix_sc_a,matrix_sc_b,matrix_overlap,para_env,context,local_data) 00097 00098 SELECT CASE (SIZE(mos)) 00099 CASE (1) 00100 s_square = 0.0_dp 00101 s_square_ideal = 0.0_dp 00102 ! let's not do this 00103 CPPrecondition(PRESENT(mo_derivs),cp_failure_level,routineP,error,failure) 00104 CASE (2) 00105 CALL get_mo_set(mo_set=mos(1)%mo_set,mo_coeff=c_alpha,homo=nalpha,uniform_occupation=uniform_occupation) 00106 CPPrecondition(uniform_occupation,cp_warning_level,routineP,error,failure) 00107 CALL get_mo_set(mo_set=mos(2)%mo_set,mo_coeff=c_beta,homo=nbeta,uniform_occupation=uniform_occupation) 00108 CPPrecondition(uniform_occupation,cp_warning_level,routineP,error,failure) 00109 CALL cp_fm_get_info(c_alpha,ncol_global=na,error=error) 00110 CALL cp_fm_get_info(c_beta,ncol_global=nb,error=error) 00111 s_square_ideal = REAL((nalpha - nbeta)*(nalpha - nbeta + 2),KIND=dp)/4.0_dp 00112 ! create overlap matrix 00113 CALL cp_fm_get_info(c_alpha,para_env=para_env,context=context,error=error) 00114 CALL cp_fm_struct_create(fm_struct_tmp,para_env=para_env,context=context, & 00115 nrow_global=na,ncol_global=nb,error=error) 00116 CALL cp_fm_create(matrix_overlap, fm_struct_tmp,name="matrix_overlap",error=error) 00117 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00118 ! create S C_beta and compute overlap 00119 CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct_tmp,nrow_global=nrow,error=error) 00120 CALL cp_fm_create(matrix_sc_b, fm_struct_tmp,name="matrix_sc_beta",error=error) 00121 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,c_beta,matrix_sc_b,nb,error=error) 00122 CALL cp_fm_gemm('T','N',na,nb,nrow,1.0_dp,c_alpha,matrix_sc_b,0.0_dp,matrix_overlap,error=error) 00123 ! invoke formula 2.271 00124 CALL cp_fm_get_info(matrix_overlap,& 00125 local_data=local_data,& 00126 nrow_local=nrow_local,& 00127 ncol_local=ncol_local,error=error) 00128 tmp=0.0_dp 00129 DO j=1,ncol_local 00130 DO i=1,nrow_local 00131 tmp=tmp+local_data(i,j)**2 00132 ENDDO 00133 ENDDO 00134 CALL mp_sum(tmp,para_env%group) 00135 s_square = s_square_ideal + nb - tmp 00136 IF (PRESENT(mo_derivs)) THEN 00137 ! this gets really wrong for fractional occupations 00138 CPPrecondition(SIZE(mo_derivs,1)==2,cp_failure_level,routineP,error,failure) 00139 CALL get_mo_set(mo_set=mos(1)%mo_set,uniform_occupation=uniform_occupation) 00140 CPPrecondition(uniform_occupation,cp_failure_level,routineP,error,failure) 00141 CALL get_mo_set(mo_set=mos(2)%mo_set,uniform_occupation=uniform_occupation) 00142 CPPrecondition(uniform_occupation,cp_failure_level,routineP,error,failure) 00143 CALL cp_fm_gemm('N','T',nrow,na,nb,-1.0_dp*strength,matrix_sc_b,matrix_overlap,1.0_dp,mo_derivs(1)%matrix,error=error) 00144 CALL cp_fm_release(matrix_sc_b,error=error) 00145 CALL cp_fm_get_info(c_alpha, matrix_struct=fm_struct_tmp,error=error) 00146 CALL cp_fm_create(matrix_sc_a, fm_struct_tmp,name="matrix_sc_alpha",error=error) 00147 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,c_alpha,matrix_sc_a,na,error=error) 00148 CALL cp_fm_gemm('N','N',nrow,nb,na,-1.0_dp*strength,matrix_sc_a,matrix_overlap,1.0_dp,mo_derivs(2)%matrix,error=error) 00149 CALL cp_fm_release(matrix_sc_a,error=error) 00150 CALL cp_fm_release(matrix_overlap,error=error) 00151 ELSE 00152 CALL cp_fm_release(matrix_sc_b,error=error) 00153 CALL cp_fm_release(matrix_overlap,error=error) 00154 ENDIF 00155 CASE DEFAULT 00156 CALL stop_program(routineN,moduleN,__LINE__,"alpha, beta, what else ?") 00157 END SELECT 00158 00159 CALL timestop(handle) 00160 00161 END SUBROUTINE compute_s_square 00162 00163 ! ***************************************************************************** 00171 SUBROUTINE s2_restraint(mos, matrix_s, mo_derivs, energy, & 00172 s2_restraint_control, just_energy, error) 00173 00174 TYPE(mo_set_p_type), DIMENSION(:), 00175 POINTER :: mos 00176 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00177 POINTER :: matrix_s 00178 TYPE(cp_fm_p_type), DIMENSION(:), 00179 POINTER :: mo_derivs 00180 REAL(kind=dp) :: energy 00181 TYPE(s2_restraint_type), POINTER :: s2_restraint_control 00182 LOGICAL :: just_energy 00183 TYPE(cp_error_type), INTENT(INOUT) :: error 00184 00185 CHARACTER(len=*), PARAMETER :: routineN = 's2_restraint', 00186 routineP = moduleN//':'//routineN 00187 00188 INTEGER :: handle 00189 LOGICAL :: failure 00190 REAL(kind=dp) :: s_square, s_square_ideal 00191 00192 CALL timeset(routineN,handle) 00193 failure = .FALSE. 00194 00195 SELECT CASE(s2_restraint_control%functional_form) 00196 CASE(do_s2_constraint) 00197 IF (just_energy) THEN 00198 CALL compute_s_square(mos, matrix_s, s_square, s_square_ideal, error=error) 00199 ELSE 00200 CALL compute_s_square(mos, matrix_s, s_square, s_square_ideal, & 00201 mo_derivs,s2_restraint_control%strength,error) 00202 ENDIF 00203 energy= s2_restraint_control%strength*(s_square-s2_restraint_control%target) 00204 s2_restraint_control%s2_order_p=s_square 00205 CASE(do_s2_restraint) ! not yet implemented 00206 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00207 CASE DEFAULT 00208 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00209 END SELECT 00210 00211 CALL timestop(handle) 00212 00213 END SUBROUTINE s2_restraint 00214 00215 END MODULE s_square_methods
1.7.3