|
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 ! ***************************************************************************** 00014 MODULE qs_resp 00015 USE atomic_charges, ONLY: print_atomic_charges 00016 USE atomic_kind_types, ONLY: atomic_kind_type,& 00017 get_atomic_kind 00018 USE cell_types, ONLY: cell_type,& 00019 pbc 00020 USE cp_output_handling, ONLY: cp_p_file,& 00021 cp_print_key_finished_output,& 00022 cp_print_key_generate_filename,& 00023 cp_print_key_should_output,& 00024 cp_print_key_unit_nr 00025 USE cp_para_types, ONLY: cp_para_env_type 00026 USE cp_subsys_types, ONLY: cp_subsys_get,& 00027 cp_subsys_type 00028 USE cp_units, ONLY: cp_unit_from_cp2k,& 00029 cp_unit_to_cp2k 00030 USE f77_blas 00031 USE input_constants, ONLY: do_resp_minus_x_dir,& 00032 do_resp_minus_y_dir,& 00033 do_resp_minus_z_dir,& 00034 do_resp_x_dir,& 00035 do_resp_y_dir,& 00036 do_resp_z_dir,& 00037 use_perd_none,& 00038 use_perd_xyz 00039 USE input_section_types, ONLY: section_get_ivals,& 00040 section_get_lval,& 00041 section_vals_get,& 00042 section_vals_get_subs_vals,& 00043 section_vals_type,& 00044 section_vals_val_get 00045 USE kinds, ONLY: default_path_length,& 00046 default_string_length,& 00047 dp 00048 USE machine, ONLY: m_flush 00049 USE mathconstants, ONLY: pi 00050 USE memory_utilities, ONLY: reallocate 00051 USE message_passing, ONLY: mp_irecv,& 00052 mp_isend,& 00053 mp_sum,& 00054 mp_wait 00055 USE particle_list_types, ONLY: particle_list_type 00056 USE particle_types, ONLY: particle_type 00057 USE pw_env_types, ONLY: pw_env_get,& 00058 pw_env_type 00059 USE pw_methods, ONLY: pw_copy,& 00060 pw_scale,& 00061 pw_transfer,& 00062 pw_zero 00063 USE pw_poisson_methods, ONLY: pw_poisson_solve 00064 USE pw_poisson_types, ONLY: pw_poisson_type 00065 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00066 pw_pool_give_back_pw,& 00067 pw_pool_type 00068 USE pw_types, ONLY: COMPLEXDATA1D,& 00069 REALDATA3D,& 00070 REALSPACE,& 00071 RECIPROCALSPACE,& 00072 pw_p_type,& 00073 pw_release,& 00074 pw_type 00075 USE qs_collocate_density, ONLY: calculate_rho_resp_all,& 00076 calculate_rho_resp_single 00077 USE qs_environment_types, ONLY: get_qs_env,& 00078 qs_environment_type 00079 USE realspace_grid_cube, ONLY: pw_to_cube 00080 USE termination, ONLY: stop_program 00081 USE timings, ONLY: timeset,& 00082 timestop 00083 #include "cp_common_uses.h" 00084 00085 IMPLICIT NONE 00086 00087 PRIVATE 00088 00089 ! *** Global parameters *** 00090 00091 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_resp' 00092 00093 PUBLIC :: resp_fit,& 00094 resp_type 00095 00096 TYPE resp_type 00097 LOGICAL :: equal_charges, itc, 00098 nonperiodic_sys, rheavies 00099 INTEGER :: nres, ncons, 00100 nrest_sec, ncons_sec, 00101 npoints, stride(3), my_fit, 00102 npoints_proc 00103 INTEGER, DIMENSION(:), POINTER :: atom_surf_list 00104 INTEGER, DIMENSION(:,:), POINTER :: fitpoints 00105 REAL(KIND=dp) :: rheavies_strength, 00106 length, eta 00107 REAL(KIND=dp), DIMENSION(3) :: box_hi, box_low 00108 REAL(KIND=dp), DIMENSION(:), POINTER :: rmin_kind, 00109 rmax_kind 00110 REAL(KIND=dp), DIMENSION(:), POINTER :: range_surf 00111 REAL(KIND=dp), DIMENSION(:), POINTER :: rhs 00112 REAL(KIND=dp), DIMENSION(:, :),POINTER :: matrix 00113 END TYPE resp_type 00114 00115 TYPE resp_p_type 00116 TYPE (resp_type), POINTER :: p_resp 00117 END TYPE resp_p_type 00118 00119 CONTAINS 00120 00121 ! ***************************************************************************** 00122 SUBROUTINE resp_fit(qs_env,error) 00123 TYPE(qs_environment_type), POINTER :: qs_env 00124 TYPE(cp_error_type), INTENT(inout) :: error 00125 00126 CHARACTER(len=*), PARAMETER :: routineN = 'resp_fit', 00127 routineP = moduleN//':'//routineN 00128 00129 INTEGER :: handle, info, my_per, natom, 00130 nvar, output_unit, stat 00131 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 00132 LOGICAL :: failure, has_resp 00133 TYPE(atomic_kind_type), DIMENSION(:), 00134 POINTER :: atomic_kind_set 00135 TYPE(cell_type), POINTER :: cell 00136 TYPE(cp_logger_type), POINTER :: logger 00137 TYPE(cp_subsys_type), POINTER :: subsys 00138 TYPE(particle_list_type), POINTER :: particles 00139 TYPE(particle_type), DIMENSION(:), 00140 POINTER :: particle_set 00141 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00142 TYPE(resp_type), POINTER :: resp_env 00143 TYPE(section_vals_type), POINTER :: cons_section, input, 00144 poisson_section, 00145 resp_section, rest_section 00146 00147 CALL timeset(routineN,handle) 00148 00149 NULLIFY(logger,atomic_kind_set,cell,subsys,particles,particle_set,input,& 00150 resp_section,cons_section,rest_section,poisson_section,resp_env,rep_sys) 00151 00152 failure=.FALSE. 00153 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00154 00155 IF (.NOT. failure) THEN 00156 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, input=input,& 00157 subsys=subsys,particle_set=particle_set, cell=cell, error=error) 00158 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 00159 error=error) 00160 CALL section_vals_get(resp_section, explicit=has_resp, error=error) 00161 END IF 00162 00163 IF (.NOT. failure .AND. has_resp) THEN 00164 logger => cp_error_get_logger(error) 00165 poisson_section => section_vals_get_subs_vals(input,"DFT%POISSON",error=error) 00166 CALL section_vals_val_get(poisson_section,"PERIODIC",i_val=my_per,error=error) 00167 CALL create_resp_type(resp_env, rep_sys, error) 00168 !initialize the RESP fitting, get all the keywords 00169 CALL init_resp(qs_env,resp_env,rep_sys,input,subsys,particle_set,atomic_kind_set,& 00170 cell,resp_section,cons_section,rest_section,error) 00171 00172 !print info 00173 CALL print_resp_parameter_info(qs_env,resp_env,rep_sys,my_per,error) 00174 00175 CALL cp_subsys_get(subsys,particles=particles,error=error) 00176 natom=particles%n_els 00177 nvar=natom+resp_env%ncons 00178 00179 ALLOCATE(ipiv(nvar),stat=stat) 00180 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00181 IF(.NOT.ASSOCIATED(resp_env%matrix)) THEN 00182 ALLOCATE(resp_env%matrix(nvar,nvar),stat=stat) 00183 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00184 ENDIF 00185 IF(.NOT.ASSOCIATED(resp_env%rhs)) THEN 00186 ALLOCATE(resp_env%rhs(nvar),stat=stat) 00187 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00188 ENDIF 00189 ipiv =0 00190 resp_env%matrix = 0.0_dp 00191 resp_env%rhs = 0.0_dp 00192 00193 ! calculate the matrix and the vector rhs 00194 SELECT CASE (my_per) 00195 CASE(use_perd_none) 00196 CALL calc_resp_matrix_nonper(qs_env,resp_env,atomic_kind_set, particles,cell,& 00197 resp_env%matrix,resp_env%rhs,natom,error) 00198 CASE(use_perd_xyz) 00199 CALL calc_resp_matrix_periodic(qs_env,resp_env,rep_sys,particles,cell,natom, error) 00200 CASE DEFAULT 00201 CALL cp_unimplemented_error(fromWhere=routineP, & 00202 message="RESP charges only implemented for nonperiodic systems"//& 00203 " or XYZ periodicity!", & 00204 error=error, error_level=cp_failure_level) 00205 END SELECT 00206 00207 output_unit=cp_print_key_unit_nr(logger,resp_section,"PRINT%PROGRAM_RUN_INFO",& 00208 extension=".resp",error=error) 00209 IF (output_unit>0) THEN 00210 WRITE(output_unit,'(T3,A,T69,I12)') "Number of potential fitting "//& 00211 "points found: ",resp_env%npoints 00212 WRITE(output_unit,'()') 00213 ENDIF 00214 00215 !adding restraints and constraints 00216 CALL add_restraints_and_constraints(qs_env,resp_env,rest_section,& 00217 subsys,natom,cons_section,particle_set,error) 00218 00219 !solve system for the values of the charges and the lagrangian multipliers 00220 CALL DGETRF(nvar,nvar,resp_env%matrix,nvar,ipiv,info) 00221 CPPrecondition(info==0,cp_failure_level,routineP,error,failure) 00222 00223 CALL DGETRS('N',nvar,1,resp_env%matrix,nvar,ipiv,resp_env%rhs,nvar,info) 00224 CPPrecondition(info==0,cp_failure_level,routineP,error,failure) 00225 00226 CALL print_resp_charges(qs_env,resp_env,output_unit,natom,error) 00227 CALL print_pot_from_resp_charges(qs_env,resp_env,particles,natom,output_unit,error) 00228 00229 DEALLOCATE(ipiv, stat=stat) 00230 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00231 CALL resp_dealloc(resp_env,rep_sys,error) 00232 CALL cp_print_key_finished_output(output_unit,logger,resp_section,& 00233 "PRINT%PROGRAM_RUN_INFO", error=error) 00234 00235 END IF 00236 00237 CALL timestop(handle) 00238 00239 END SUBROUTINE resp_fit 00240 !**************************************************************************** 00241 !\brief creates the resp_type structure 00242 !\param error variable to control error logging, stopping,... 00243 ! see module cp_error_handling 00244 !**************************************************************************** 00245 SUBROUTINE create_resp_type(resp_env, rep_sys, error) 00246 TYPE(resp_type), POINTER :: resp_env 00247 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00248 TYPE(cp_error_type), INTENT(inout) :: error 00249 00250 CHARACTER(len=*), PARAMETER :: routineN = 'create_resp_type', 00251 routineP = moduleN//':'//routineN 00252 00253 INTEGER :: stat 00254 LOGICAL :: failure 00255 00256 failure = .FALSE. 00257 IF (ASSOCIATED(resp_env)) CALL resp_dealloc(resp_env,rep_sys,error) 00258 ALLOCATE(resp_env, stat=stat) 00259 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00260 00261 IF (.NOT. failure) THEN 00262 NULLIFY(resp_env%matrix,& 00263 resp_env%fitpoints,& 00264 resp_env%rmin_kind,& 00265 resp_env%rmax_kind,& 00266 resp_env%rhs) 00267 00268 ENDIF 00269 00270 resp_env%equal_charges=.FALSE. 00271 resp_env%itc=.FALSE. 00272 resp_env%nonperiodic_sys=.FALSE. 00273 resp_env%rheavies=.FALSE. 00274 00275 resp_env%box_hi=0.0_dp 00276 resp_env%box_low=0.0_dp 00277 00278 resp_env%ncons=0 00279 resp_env%ncons_sec=0 00280 resp_env%nres=0 00281 resp_env%nrest_sec=0 00282 resp_env%npoints=0 00283 resp_env%npoints_proc=0 00284 00285 END SUBROUTINE create_resp_type 00286 !**************************************************************************** 00287 !\brief deallocates the resp_type structure 00288 !\param error variable to control error logging, stopping,... 00289 ! see module cp_error_handling 00290 !**************************************************************************** 00291 SUBROUTINE resp_dealloc(resp_env,rep_sys, error) 00292 TYPE(resp_type), POINTER :: resp_env 00293 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00294 TYPE(cp_error_type), INTENT(inout) :: error 00295 00296 CHARACTER(len=*), PARAMETER :: routineN = 'resp_dealloc', 00297 routineP = moduleN//':'//routineN 00298 00299 INTEGER :: i, stat 00300 00301 IF (ASSOCIATED(resp_env)) THEN 00302 IF (ASSOCIATED(resp_env%matrix)) THEN 00303 DEALLOCATE(resp_env%matrix, stat=stat) 00304 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00305 ENDIF 00306 IF (ASSOCIATED(resp_env%rhs)) THEN 00307 DEALLOCATE(resp_env%rhs, stat=stat) 00308 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00309 ENDIF 00310 IF (ASSOCIATED(resp_env%fitpoints)) THEN 00311 DEALLOCATE(resp_env%fitpoints, stat=stat) 00312 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00313 ENDIF 00314 IF (ASSOCIATED(resp_env%rmin_kind))THEN 00315 DEALLOCATE(resp_env%rmin_kind, stat=stat) 00316 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00317 ENDIF 00318 IF (ASSOCIATED(resp_env%rmax_kind))THEN 00319 DEALLOCATE(resp_env%rmax_kind, stat=stat) 00320 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00321 ENDIF 00322 DEALLOCATE(resp_env, stat=stat) 00323 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00324 ENDIF 00325 IF (ASSOCIATED(rep_sys)) THEN 00326 DO i=1, SIZE(rep_sys) 00327 DEALLOCATE(rep_sys(i)%p_resp%atom_surf_list, stat=stat) 00328 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00329 DEALLOCATE(rep_sys(i)%p_resp, stat=stat) 00330 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00331 ENDDO 00332 DEALLOCATE(rep_sys, stat=stat) 00333 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00334 ENDIF 00335 00336 END SUBROUTINE resp_dealloc 00337 !**************************************************************************** 00338 !\brief intializes the resp fit. Getting the parameters 00339 !**************************************************************************** 00340 SUBROUTINE init_resp(qs_env,resp_env,rep_sys,input,subsys,particle_set,& 00341 atomic_kind_set,cell,resp_section,cons_section,rest_section,error) 00342 00343 TYPE(qs_environment_type), POINTER :: qs_env 00344 TYPE(resp_type), POINTER :: resp_env 00345 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00346 TYPE(section_vals_type), POINTER :: input 00347 TYPE(cp_subsys_type), POINTER :: subsys 00348 TYPE(particle_type), DIMENSION(:), 00349 POINTER :: particle_set 00350 TYPE(atomic_kind_type), DIMENSION(:), 00351 POINTER :: atomic_kind_set 00352 TYPE(cell_type), POINTER :: cell 00353 TYPE(section_vals_type), POINTER :: resp_section, cons_section, 00354 rest_section 00355 TYPE(cp_error_type), INTENT(inout) :: error 00356 00357 CHARACTER(len=*), PARAMETER :: routineN = 'init_resp', 00358 routineP = moduleN//':'//routineN 00359 00360 INTEGER :: handle, i, nrep, stat 00361 INTEGER, DIMENSION(:), POINTER :: atom_list_cons, my_stride 00362 LOGICAL :: explicit, failure 00363 TYPE(section_vals_type), POINTER :: nonperiodic_section, 00364 periodic_section 00365 00366 CALL timeset(routineN,handle) 00367 00368 NULLIFY(atom_list_cons, my_stride, nonperiodic_section, periodic_section) 00369 failure=.FALSE. 00370 00371 ! get the subsections 00372 nonperiodic_section=>section_vals_get_subs_vals(resp_section,"NONPERIODIC_SYS",& 00373 error=error) 00374 periodic_section=>section_vals_get_subs_vals(resp_section,"PERIODIC_SYS",& 00375 error=error) 00376 cons_section=>section_vals_get_subs_vals(resp_section,"CONSTRAINT",& 00377 error=error) 00378 rest_section=>section_vals_get_subs_vals(resp_section,"RESTRAINT",& 00379 error=error) 00380 00381 ! get and set the parameters for nonperiodic (non-surface) systems 00382 CALL get_parameter_nonperiodic_sys(resp_env, nonperiodic_section, cell, & 00383 atomic_kind_set, error) 00384 00385 ! get the parameter for periodic/surface systems 00386 CALL section_vals_get(periodic_section, explicit=explicit, n_repetition=nrep, error=error) 00387 IF(explicit) THEN 00388 ALLOCATE(rep_sys(nrep), stat=stat) 00389 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00390 DO i=1, nrep 00391 ALLOCATE(rep_sys(i)%p_resp, stat=stat) 00392 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00393 NULLIFY(rep_sys(i)%p_resp%range_surf, rep_sys(i)%p_resp%atom_surf_list) 00394 CALL section_vals_val_get(periodic_section,"RANGE",r_vals=rep_sys(i)%p_resp%range_surf,& 00395 i_rep_section=i, error=error) 00396 CALL section_vals_val_get(periodic_section,"LENGTH", r_val=rep_sys(i)%p_resp%length,& 00397 i_rep_section=i, error=error) 00398 CALL section_vals_val_get(periodic_section,"SURF_DIRECTION",& 00399 i_rep_section=i, i_val=rep_sys(i)%p_resp%my_fit,error=error) 00400 IF(ANY(rep_sys(i)%p_resp%range_surf<0.0_dp)) THEN 00401 CALL stop_program(routineN,moduleN,__LINE__,& 00402 "Numbers in RANGE in PERIODIC_SYS cannot be negative.") 00403 ENDIF 00404 IF(rep_sys(i)%p_resp%length<=EPSILON(0.0_dp)) THEN 00405 CALL stop_program(routineN,moduleN,__LINE__,& 00406 "Parameter LENGTH in PERIODIC_SYS has to be larger than zero.") 00407 ENDIF 00408 !list of atoms specifing the surface 00409 CALL build_atom_list(periodic_section,subsys,rep_sys(i)%p_resp%atom_surf_list,rep=i,error=error) 00410 ENDDO 00411 ENDIF 00412 00413 00414 ! get the parameters for the constraint and restraint sections 00415 CALL section_vals_get(cons_section, explicit=explicit, error=error) 00416 IF (explicit) THEN 00417 CALL section_vals_get(cons_section,n_repetition=resp_env%ncons_sec,error=error) 00418 DO i=1,resp_env%ncons_sec 00419 CALL section_vals_val_get(cons_section,"EQUAL_CHARGES",l_val=resp_env%equal_charges,& 00420 explicit=explicit,error=error) 00421 IF(.NOT.explicit) CYCLE 00422 CALL build_atom_list(cons_section,subsys,atom_list_cons,i,error=error) 00423 !instead of using EQUAL_CHARGES the constraint sections could be repeated 00424 resp_env%ncons=resp_env%ncons+SIZE(atom_list_cons)-2 00425 DEALLOCATE(atom_list_cons,stat=stat) 00426 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00427 ENDDO 00428 ENDIF 00429 CALL section_vals_get(rest_section, explicit=explicit, error=error) 00430 IF (explicit) THEN 00431 CALL section_vals_get(rest_section,n_repetition=resp_env%nrest_sec,error=error) 00432 ENDIF 00433 resp_env%ncons=resp_env%ncons+resp_env%ncons_sec 00434 resp_env%nres=resp_env%nres+resp_env%nrest_sec 00435 00436 ! get the general keywords 00437 CALL section_vals_val_get(resp_section,"INTEGER_TOTAL_CHARGE",& 00438 l_val=resp_env%itc,error=error) 00439 IF (resp_env%itc) resp_env%ncons=resp_env%ncons+1 00440 00441 CALL section_vals_val_get(resp_section,"RESTRAIN_HEAVIES_TO_ZERO",& 00442 l_val=resp_env%rheavies,error=error) 00443 IF (resp_env%rheavies) THEN 00444 CALL section_vals_val_get(resp_section,"RESTRAIN_HEAVIES_STRENGTH",& 00445 r_val=resp_env%rheavies_strength,error=error) 00446 ENDIF 00447 CALL section_vals_val_get(resp_section,"STRIDE",i_vals=my_stride,error=error) 00448 CALL cp_assert(SIZE(my_stride)==1.OR.SIZE(my_stride)==3,cp_fatal_level,cp_assertion_failed,routineP,& 00449 "STRIDE keyword can accept only 1 (the same for X,Y,Z) or 3 values. Correct your input file."//& 00450 CPSourceFileRef,& 00451 only_ionode=.TRUE.) 00452 IF (SIZE(my_stride)==1) THEN 00453 DO i = 1,3 00454 resp_env%stride(i) = my_stride(1) 00455 END DO 00456 ELSE 00457 resp_env%stride = my_stride(1:3) 00458 END IF 00459 CALL section_vals_val_get(resp_section,"WIDTH", r_val=resp_env%eta, error=error) 00460 00461 CALL timestop(handle) 00462 00463 END SUBROUTINE init_resp 00464 !**************************************************************************** 00465 !\brief getting the parameters for nonperiodic/non-surface systems 00466 !**************************************************************************** 00467 SUBROUTINE get_parameter_nonperiodic_sys(resp_env,nonperiodic_section,cell,& 00468 atomic_kind_set,error) 00469 00470 TYPE(resp_type), POINTER :: resp_env 00471 TYPE(section_vals_type), POINTER :: nonperiodic_section 00472 TYPE(cell_type), POINTER :: cell 00473 TYPE(atomic_kind_type), DIMENSION(:), 00474 POINTER :: atomic_kind_set 00475 TYPE(cp_error_type), INTENT(inout) :: error 00476 00477 CHARACTER(len=*), PARAMETER :: 00478 routineN = 'get_parameter_nonperiodic_sys', 00479 routineP = moduleN//':'//routineN 00480 00481 CHARACTER(LEN=2) :: symbol 00482 CHARACTER(LEN=default_string_length), 00483 DIMENSION(:), POINTER :: tmpstringlist 00484 INTEGER :: ikind, j, kind_number, nkind, 00485 nrep_rmax, nrep_rmin, stat 00486 LOGICAL :: explicit, failure 00487 REAL(KIND=dp) :: rmax, rmin 00488 TYPE(atomic_kind_type), POINTER :: atomic_kind 00489 00490 failure=.FALSE. 00491 nrep_rmin=0 00492 nrep_rmax=0 00493 nkind=SIZE(atomic_kind_set) 00494 00495 CALL section_vals_get(nonperiodic_section, explicit=explicit, error=error) 00496 IF(explicit) THEN 00497 resp_env%nonperiodic_sys=.TRUE. 00498 CALL section_vals_val_get(nonperiodic_section,"RMIN",r_val=rmin,& 00499 error=error) 00500 CALL section_vals_val_get(nonperiodic_section,"RMAX",r_val=rmax,& 00501 error=error) 00502 CALL section_vals_val_get(nonperiodic_section,"RMIN_KIND", n_rep_val=nrep_rmin, error=error) 00503 CALL section_vals_val_get(nonperiodic_section,"RMAX_KIND", n_rep_val=nrep_rmax, error=error) 00504 ALLOCATE(resp_env%rmin_kind(nkind),stat=stat) 00505 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00506 ALLOCATE(resp_env%rmax_kind(nkind),stat=stat) 00507 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00508 resp_env%rmin_kind(:)=rmin 00509 resp_env%rmax_kind(:)=rmax 00510 DO j=1, nrep_rmin 00511 CALL section_vals_val_get(nonperiodic_section,"RMIN_KIND", i_rep_val=j,& 00512 c_vals=tmpstringlist, error=error) 00513 DO ikind=1,nkind 00514 atomic_kind => atomic_kind_set(ikind) 00515 CALL get_atomic_kind(atomic_kind,element_symbol=symbol,kind_number=kind_number) 00516 IF(TRIM(tmpstringlist(2))==TRIM(symbol)) THEN 00517 READ(tmpstringlist(1),*) resp_env%rmin_kind(kind_number) 00518 resp_env%rmin_kind(kind_number)=cp_unit_to_cp2k(resp_env%rmin_kind(kind_number),& 00519 "angstrom",error=error) 00520 ENDIF 00521 ENDDO 00522 ENDDO 00523 DO j=1, nrep_rmax 00524 CALL section_vals_val_get(nonperiodic_section,"RMAX_KIND", i_rep_val=j,& 00525 c_vals=tmpstringlist, error=error) 00526 DO ikind=1,nkind 00527 atomic_kind => atomic_kind_set(ikind) 00528 CALL get_atomic_kind(atomic_kind,element_symbol=symbol,kind_number=kind_number) 00529 IF(TRIM(tmpstringlist(2))==TRIM(symbol)) THEN 00530 READ(tmpstringlist(1),*) resp_env%rmax_kind(kind_number) 00531 resp_env%rmax_kind(kind_number)=cp_unit_to_cp2k(resp_env%rmax_kind(kind_number),& 00532 "angstrom",error=error) 00533 ENDIF 00534 ENDDO 00535 ENDDO 00536 00537 resp_env%box_hi=(/cell%hmat(1,1),cell%hmat(2,2),cell%hmat(3,3)/) 00538 resp_env%box_low=0.0_dp 00539 CALL section_vals_val_get(nonperiodic_section,"X_HI",explicit=explicit,error=error) 00540 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"X_HI",& 00541 r_val=resp_env%box_hi(1),error=error) 00542 CALL section_vals_val_get(nonperiodic_section,"X_LOW",explicit=explicit,error=error) 00543 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"X_LOW",& 00544 r_val=resp_env%box_low(1),error=error) 00545 CALL section_vals_val_get(nonperiodic_section,"Y_HI",explicit=explicit,error=error) 00546 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"Y_HI",& 00547 r_val=resp_env%box_hi(2),error=error) 00548 CALL section_vals_val_get(nonperiodic_section,"Y_LOW",explicit=explicit,error=error) 00549 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"Y_LOW",& 00550 r_val=resp_env%box_low(2),error=error) 00551 CALL section_vals_val_get(nonperiodic_section,"Z_HI",explicit=explicit,error=error) 00552 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"Z_HI",& 00553 r_val=resp_env%box_hi(3),error=error) 00554 CALL section_vals_val_get(nonperiodic_section,"Z_LOW",explicit=explicit,error=error) 00555 IF (explicit) CALL section_vals_val_get(nonperiodic_section,"Z_LOW",& 00556 r_val=resp_env%box_low(3),error=error) 00557 ENDIF 00558 00559 END SUBROUTINE get_parameter_nonperiodic_sys 00560 !**************************************************************************** 00561 !\brief building atom lists for different sections of RESP 00562 !**************************************************************************** 00563 SUBROUTINE build_atom_list(section,subsys,atom_list,rep,error) 00564 00565 TYPE(section_vals_type), POINTER :: section 00566 TYPE(cp_subsys_type), POINTER :: subsys 00567 INTEGER, DIMENSION(:), POINTER :: atom_list 00568 INTEGER, INTENT(IN), OPTIONAL :: rep 00569 TYPE(cp_error_type), INTENT(inout) :: error 00570 00571 CHARACTER(len=*), PARAMETER :: routineN = 'build_atom_list', 00572 routineP = moduleN//':'//routineN 00573 00574 INTEGER :: atom_a, atom_b, handle, i, 00575 irep, j, max_index, n_var, 00576 num_atom, stat 00577 INTEGER, DIMENSION(:), POINTER :: indexes 00578 LOGICAL :: failure, index_in_range 00579 00580 CALL timeset(routineN,handle) 00581 00582 NULLIFY(indexes) 00583 failure=.FALSE. 00584 irep=1 00585 IF(PRESENT(rep)) irep=rep 00586 00587 CALL section_vals_val_get(section,"ATOM_LIST",i_rep_section=irep,& 00588 n_rep_val=n_var,error=error) 00589 num_atom=0 00590 DO i=1,n_var 00591 CALL section_vals_val_get(section,"ATOM_LIST",i_rep_section=irep,& 00592 i_rep_val=i,i_vals=indexes,error=error) 00593 num_atom=num_atom + SIZE(indexes) 00594 ENDDO 00595 ALLOCATE(atom_list(num_atom),stat=stat) 00596 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00597 atom_list=0 00598 num_atom=1 00599 DO i=1,n_var 00600 CALL section_vals_val_get(section,"ATOM_LIST",i_rep_section=irep,& 00601 i_rep_val=i,i_vals=indexes,error=error) 00602 atom_list(num_atom:num_atom+SIZE(indexes)-1)=indexes(:) 00603 num_atom = num_atom + SIZE(indexes) 00604 ENDDO 00605 !check atom list 00606 num_atom=num_atom-1 00607 max_index=SIZE(subsys%particles%els) 00608 CPPrecondition(SIZE(atom_list) /= 0,cp_failure_level,routineP,error,failure) 00609 index_in_range=(MAXVAL(atom_list)<= max_index)& 00610 .AND.(MINVAL(atom_list) > 0) 00611 CPPostcondition(index_in_range,cp_failure_level,routineP,error,failure) 00612 DO i=1,num_atom 00613 DO j=i+1,num_atom 00614 atom_a=atom_list(i) 00615 atom_b=atom_list(j) 00616 IF(atom_a==atom_b) & 00617 CALL stop_program(routineN,moduleN,__LINE__,& 00618 "There are atoms doubled in atom list for RESP.") 00619 ENDDO 00620 ENDDO 00621 00622 CALL timestop(handle) 00623 00624 END SUBROUTINE build_atom_list 00625 00626 !**************************************************************************** 00627 !\brief build matrix and vector for nonperiodic RESP fitting 00628 !**************************************************************************** 00629 SUBROUTINE calc_resp_matrix_nonper(qs_env,resp_env,atomic_kind_set,particles,& 00630 cell,matrix,rhs,natom,error) 00631 00632 TYPE(qs_environment_type), POINTER :: qs_env 00633 TYPE(resp_type), POINTER :: resp_env 00634 TYPE(atomic_kind_type), DIMENSION(:), 00635 POINTER :: atomic_kind_set 00636 TYPE(particle_list_type), POINTER :: particles 00637 TYPE(cell_type), POINTER :: cell 00638 REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix 00639 REAL(KIND=dp), DIMENSION(:), POINTER :: rhs 00640 INTEGER, INTENT(IN) :: natom 00641 TYPE(cp_error_type), INTENT(inout) :: error 00642 00643 CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_matrix_nonper', 00644 routineP = moduleN//':'//routineN 00645 00646 INTEGER :: bo(2,3), gbo(2,3), handle, i, 00647 ikind, jx, jy, jz, k, 00648 kind_number, l, m, nkind, 00649 now, np(3), p, stat 00650 LOGICAL :: failure 00651 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range 00652 REAL(KIND=dp) :: dh(3,3), dvol, r(3), rmax, 00653 rmin, vec(3), vec_pbc(3), vj 00654 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dist 00655 TYPE(cp_logger_type), POINTER :: logger 00656 TYPE(particle_type), DIMENSION(:), 00657 POINTER :: particle_set 00658 TYPE(pw_type), POINTER :: v_hartree_pw 00659 TYPE(section_vals_type), POINTER :: input, print_key, resp_section 00660 00661 CALL timeset(routineN,handle) 00662 00663 NULLIFY(logger,input,particle_set,print_key,resp_section,v_hartree_pw) 00664 failure=.FALSE. 00665 00666 IF (.NOT.cell%orthorhombic) THEN 00667 CALL cp_unimplemented_error(fromWhere=routineP, & 00668 message="Nonperiodic solution for RESP charges only"//& 00669 " implemented for orthorhombic cells!", & 00670 error=error, error_level=cp_failure_level) 00671 END IF 00672 IF(.NOT.resp_env%nonperiodic_sys) THEN 00673 CALL stop_program(routineN,moduleN,__LINE__,& 00674 "Nonperiodic solution for RESP charges (i.e. nonperiodic"//& 00675 " Poisson solver) can only be used with section NONPERIODC_SYS") 00676 ENDIF 00677 CALL get_qs_env(qs_env=qs_env,input=input,particle_set=particle_set,error=error) 00678 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 00679 error=error) 00680 print_key => section_vals_get_subs_vals(resp_section,& 00681 "PRINT%V_RESP_CUBE",error=error) 00682 logger => cp_error_get_logger(error) 00683 00684 v_hartree_pw => qs_env%ks_env%v_hartree_rspace%pw 00685 bo=v_hartree_pw%pw_grid%bounds_local 00686 gbo=v_hartree_pw%pw_grid%bounds 00687 np=v_hartree_pw%pw_grid%npts 00688 dh=v_hartree_pw%pw_grid%dh 00689 dvol=v_hartree_pw%pw_grid%dvol 00690 nkind=SIZE(atomic_kind_set) 00691 00692 ALLOCATE(dist(natom),stat=stat) 00693 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00694 ALLOCATE(not_in_range(natom,2),stat=stat) 00695 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00696 00697 !only store fitting points when printing RESP potential 00698 IF (BTEST(cp_print_key_should_output(logger%iter_info,resp_section,& 00699 "PRINT%V_RESP_CUBE",error=error),cp_p_file)) THEN 00700 IF(.NOT.ASSOCIATED(resp_env%fitpoints)) THEN 00701 now = 1000 00702 ALLOCATE(resp_env%fitpoints(3,now),stat=stat) 00703 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00704 ELSE 00705 now = SIZE(resp_env%fitpoints,2) 00706 END IF 00707 ENDIF 00708 00709 DO jz=bo(1,3),bo(2,3) 00710 DO jy=bo(1,2),bo(2,2) 00711 DO jx=bo(1,1),bo(2,1) 00712 IF (.NOT.(MODULO(jz,resp_env%stride(3))==0)) CYCLE 00713 IF (.NOT.(MODULO(jy,resp_env%stride(2))==0)) CYCLE 00714 IF (.NOT.(MODULO(jx,resp_env%stride(1))==0)) CYCLE 00715 !bounds bo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0) 00716 l=jx - gbo(1,1) 00717 k=jy - gbo(1,2) 00718 p=jz - gbo(1,3) 00719 r(3)=p*dh(3,3)+k*dh(3,2)+l*dh(3,1) 00720 r(2)=p*dh(2,3)+k*dh(2,2)+l*dh(2,1) 00721 r(1)=p*dh(1,3)+k*dh(1,2)+l*dh(1,1) 00722 IF (r(3)<resp_env%box_low(3).OR.r(3)>resp_env%box_hi(3)) CYCLE 00723 IF (r(2)<resp_env%box_low(2).OR.r(2)>resp_env%box_hi(2)) CYCLE 00724 IF (r(1)<resp_env%box_low(1).OR.r(1)>resp_env%box_hi(1)) CYCLE 00725 ! compute distance from the grid point to all atoms 00726 not_in_range=.FALSE. 00727 DO i=1,natom 00728 vec=r-particles%els(i)%r 00729 vec_pbc(1) = vec(1) - cell%hmat(1,1)*ANINT(cell%h_inv(1,1)*vec(1)) 00730 vec_pbc(2) = vec(2) - cell%hmat(2,2)*ANINT(cell%h_inv(2,2)*vec(2)) 00731 vec_pbc(3) = vec(3) - cell%hmat(3,3)*ANINT(cell%h_inv(3,3)*vec(3)) 00732 dist(i)=SQRT(SUM(vec_pbc**2)) 00733 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind,& 00734 kind_number=kind_number) 00735 DO ikind=1,nkind 00736 IF(ikind==kind_number) THEN 00737 rmin=resp_env%rmin_kind(ikind) 00738 rmax=resp_env%rmax_kind(ikind) 00739 EXIT 00740 ENDIF 00741 ENDDO 00742 IF(dist(i)<rmin+1.0E-13_dp) not_in_range(i,1)=.TRUE. 00743 IF(dist(i)>rmax-1.0E-13_dp) not_in_range(i,2)=.TRUE. 00744 ENDDO 00745 ! check if the point is sufficiently close and far. if OK, we can use 00746 ! the point for fitting, add/subtract 1.0E-13 to get rid of rounding errors when shifting atoms 00747 IF(ANY(not_in_range(:,1)).OR.ALL(not_in_range(:,2))) CYCLE 00748 resp_env%npoints = resp_env%npoints + 1 00749 IF (BTEST(cp_print_key_should_output(logger%iter_info,resp_section,& 00750 "PRINT%V_RESP_CUBE",error=error),cp_p_file)) THEN 00751 IF(resp_env%npoints > now) THEN 00752 now = 2*now 00753 CALL reallocate(resp_env%fitpoints,1,3,1,now) 00754 ENDIF 00755 resp_env%fitpoints(1,resp_env%npoints) = jx 00756 resp_env%fitpoints(2,resp_env%npoints) = jy 00757 resp_env%fitpoints(3,resp_env%npoints) = jz 00758 ENDIF 00759 ! correct for the fact that v_hartree is scaled by dvol, and has the opposite sign 00760 IF (qs_env%qmmm) THEN 00761 ! If it's a QM/MM run let's remove the contribution of the MM potential out of the Hartree pot 00762 vj=-v_hartree_pw%cr3d(jx,jy,jz)/dvol+qs_env%ks_qmmm_env%v_qmmm_rspace%pw%cr3d(jx,jy,jz) 00763 ELSE 00764 vj=-v_hartree_pw%cr3d(jx,jy,jz)/dvol 00765 END IF 00766 dist=1.0_dp/dist 00767 00768 DO i=1,natom 00769 DO m=1,natom 00770 matrix(m,i)=matrix(m,i)+2.0_dp*dist(i)*dist(m) 00771 ENDDO 00772 rhs(i)=rhs(i)+2.0_dp*vj*dist(i) 00773 ENDDO 00774 ENDDO 00775 ENDDO 00776 ENDDO 00777 00778 resp_env%npoints_proc=resp_env%npoints 00779 CALL mp_sum(resp_env%npoints,v_hartree_pw%pw_grid%para%group) 00780 CALL mp_sum(matrix,v_hartree_pw%pw_grid%para%group) 00781 CALL mp_sum(rhs,v_hartree_pw%pw_grid%para%group) 00782 !weighted sum 00783 matrix=matrix/resp_env%npoints 00784 rhs=rhs/resp_env%npoints 00785 00786 DEALLOCATE(dist,stat=stat) 00787 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00788 DEALLOCATE(not_in_range,stat=stat) 00789 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00790 00791 CALL timestop(handle) 00792 00793 END SUBROUTINE calc_resp_matrix_nonper 00794 !**************************************************************************** 00795 !\brief build matrix and vector for periodic RESP fitting 00796 !**************************************************************************** 00797 SUBROUTINE calc_resp_matrix_periodic(qs_env,resp_env,rep_sys,particles,cell,& 00798 natom,error) 00799 00800 TYPE(qs_environment_type), POINTER :: qs_env 00801 TYPE(resp_type), POINTER :: resp_env 00802 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00803 TYPE(particle_list_type), POINTER :: particles 00804 TYPE(cell_type), POINTER :: cell 00805 INTEGER, INTENT(IN) :: natom 00806 TYPE(cp_error_type), INTENT(inout) :: error 00807 00808 CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_matrix_periodic', 00809 routineP = moduleN//':'//routineN 00810 00811 INTEGER :: handle, i, ip, j, jx, jy, jz, 00812 stat 00813 LOGICAL :: failure 00814 REAL(KIND=dp) :: normalize_factor 00815 REAL(KIND=dp), ALLOCATABLE, 00816 DIMENSION(:, :) :: vpot 00817 TYPE(cp_para_env_type), POINTER :: para_env 00818 TYPE(pw_env_type), POINTER :: pw_env 00819 TYPE(pw_p_type) :: rho_ga, va_gspace, va_rspace 00820 TYPE(pw_poisson_type), POINTER :: poisson_env 00821 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00822 00823 CALL timeset(routineN,handle) 00824 00825 failure=.FALSE. 00826 NULLIFY(pw_env, para_env,auxbas_pw_pool,poisson_env) 00827 00828 IF(.NOT.ALL(cell%perd/=0)) THEN 00829 CALL stop_program(routineN,moduleN,__LINE__,& 00830 "Periodic solution for RESP (with periodic Poisson solver)"//& 00831 " can only be obtained with a cell that has XYZ periodicity") 00832 ENDIF 00833 00834 CALL get_qs_env(qs_env, pw_env=pw_env,para_env=para_env,& 00835 error=error) 00836 00837 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 00838 poisson_env=poisson_env, error=error) 00839 CALL pw_pool_create_pw(auxbas_pw_pool,& 00840 rho_ga%pw,& 00841 use_data=COMPLEXDATA1D,& 00842 in_space=RECIPROCALSPACE,& 00843 error=error) 00844 CALL pw_pool_create_pw(auxbas_pw_pool,& 00845 va_gspace%pw,& 00846 use_data=COMPLEXDATA1D,& 00847 in_space=RECIPROCALSPACE,& 00848 error=error) 00849 CALL pw_pool_create_pw(auxbas_pw_pool,& 00850 va_rspace%pw,& 00851 use_data=REALDATA3D,& 00852 in_space=REALSPACE,& 00853 error=error) 00854 00855 !get fitting points and store them in resp_env%fitpoints 00856 CALL get_fitting_points(qs_env,resp_env,rep_sys,particles=particles,& 00857 cell=cell,error=error) 00858 ALLOCATE(vpot(resp_env%npoints,natom),stat=stat) 00859 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00860 normalize_factor=SQRT((resp_env%eta/pi)**3) 00861 00862 DO i=1,natom 00863 !collocate gaussian for each atom 00864 CALL pw_zero(rho_ga%pw, error=error) 00865 CALL calculate_rho_resp_single(rho_ga,qs_env,resp_env%eta,i,error) 00866 !calculate potential va and store the part needed for fitting in vpot 00867 CALL pw_zero(va_gspace%pw, error=error) 00868 CALL pw_poisson_solve(poisson_env,rho_ga%pw,vhartree=va_gspace%pw,error=error) 00869 CALL pw_zero(va_rspace%pw, error=error) 00870 CALL pw_transfer(va_gspace%pw,va_rspace%pw,error=error) 00871 CALL pw_scale(va_rspace%pw,normalize_factor,error=error) 00872 DO ip=1,resp_env%npoints 00873 jx = resp_env%fitpoints(1,ip) 00874 jy = resp_env%fitpoints(2,ip) 00875 jz = resp_env%fitpoints(3,ip) 00876 vpot(ip,i) = va_rspace%pw%cr3d(jx,jy,jz) 00877 END DO 00878 ENDDO 00879 00880 CALL pw_release(va_gspace%pw,error=error) 00881 CALL pw_release(va_rspace%pw,error=error) 00882 CALL pw_release(rho_ga%pw,error=error) 00883 00884 DO i=1,natom 00885 DO j=1,natom 00886 ! calculate matrix 00887 resp_env%matrix(i,j)=resp_env%matrix(i,j) + 2.0_dp*SUM(vpot(:,i)*vpot(:,j)) 00888 ENDDO 00889 ! calculate vector resp_env%rhs 00890 CALL calculate_rhs(qs_env,resp_env,resp_env%rhs(i),vpot(:,i),error=error) 00891 ENDDO 00892 00893 DEALLOCATE(vpot,stat=stat) 00894 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 00895 00896 resp_env%npoints_proc=resp_env%npoints 00897 CALL mp_sum(resp_env%npoints,para_env%group) 00898 CALL mp_sum(resp_env%matrix,para_env%group) 00899 CALL mp_sum(resp_env%rhs,para_env%group) 00900 !weighted sum 00901 resp_env%matrix=resp_env%matrix/resp_env%npoints 00902 resp_env%rhs=resp_env%rhs/resp_env%npoints 00903 00904 CALL timestop(handle) 00905 00906 END SUBROUTINE calc_resp_matrix_periodic 00907 00908 !**************************************************************************** 00909 !\brief get RESP fitting points for the periodic fitting 00910 !**************************************************************************** 00911 SUBROUTINE get_fitting_points(qs_env,resp_env,rep_sys,particles,cell,error) 00912 00913 TYPE(qs_environment_type), POINTER :: qs_env 00914 TYPE(resp_type), POINTER :: resp_env 00915 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 00916 TYPE(particle_list_type), POINTER :: particles 00917 TYPE(cell_type), POINTER :: cell 00918 TYPE(cp_error_type), INTENT(inout) :: error 00919 00920 CHARACTER(len=*), PARAMETER :: routineN = 'get_fitting_points', 00921 routineP = moduleN//':'//routineN 00922 00923 INTEGER :: bo(2,3), gbo(2,3), handle, i, iatom, ikind, in_x, in_y, in_z, 00924 jx, jy, jz, k, kind_number, l, m, natom, nkind, now, output_unit, p, 00925 stat 00926 LOGICAL :: failure 00927 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range 00928 REAL(KIND=dp) :: dh(3,3), r(3), rmax, rmin, 00929 vec_pbc(3) 00930 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dist 00931 TYPE(atomic_kind_type), DIMENSION(:), 00932 POINTER :: atomic_kind_set 00933 TYPE(cp_logger_type), POINTER :: logger 00934 TYPE(cp_para_env_type), POINTER :: para_env 00935 TYPE(particle_type), DIMENSION(:), 00936 POINTER :: particle_set 00937 TYPE(pw_type), POINTER :: v_hartree_pw 00938 TYPE(section_vals_type), POINTER :: input, resp_section 00939 00940 CALL timeset(routineN,handle) 00941 00942 failure=.FALSE. 00943 NULLIFY(atomic_kind_set,v_hartree_pw,para_env,logger,input,particle_set,& 00944 resp_section) 00945 00946 CALL get_qs_env(qs_env,input=input,particle_set=particle_set,& 00947 atomic_kind_set=atomic_kind_set,para_env=para_env,& 00948 error=error) 00949 00950 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 00951 error=error) 00952 logger => cp_error_get_logger(error) 00953 output_unit=cp_print_key_unit_nr(logger,resp_section,& 00954 "PRINT%COORD_FIT_POINTS",& 00955 extension=".xyz",& 00956 file_status="REPLACE",& 00957 file_action="WRITE",& 00958 file_form="FORMATTED",& 00959 error=error) 00960 00961 v_hartree_pw => qs_env%ks_env%v_hartree_rspace%pw 00962 bo=v_hartree_pw%pw_grid%bounds_local 00963 gbo=v_hartree_pw%pw_grid%bounds 00964 dh=v_hartree_pw%pw_grid%dh 00965 natom=SIZE(particles%els) 00966 nkind=SIZE(atomic_kind_set) 00967 00968 IF(.NOT.ASSOCIATED(resp_env%fitpoints)) THEN 00969 now = 1000 00970 ALLOCATE(resp_env%fitpoints(3,now),stat=stat) 00971 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00972 ELSE 00973 now = SIZE(resp_env%fitpoints,2) 00974 END IF 00975 00976 ALLOCATE(dist(natom),stat=stat) 00977 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00978 ALLOCATE(not_in_range(natom,2),stat=stat) 00979 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 00980 00981 !every proc gets another bo, grid is distributed 00982 DO jz=bo(1,3),bo(2,3) 00983 IF (.NOT.(MODULO(jz,resp_env%stride(3))==0)) CYCLE 00984 DO jy=bo(1,2),bo(2,2) 00985 IF (.NOT.(MODULO(jy,resp_env%stride(2))==0)) CYCLE 00986 DO jx=bo(1,1),bo(2,1) 00987 IF (.NOT.(MODULO(jx,resp_env%stride(1))==0)) CYCLE 00988 !bounds gbo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0) 00989 l=jx - gbo(1,1) 00990 k=jy - gbo(1,2) 00991 p=jz - gbo(1,3) 00992 r(3)=p*dh(3,3)+k*dh(3,2)+l*dh(3,1) 00993 r(2)=p*dh(2,3)+k*dh(2,2)+l*dh(2,1) 00994 r(1)=p*dh(1,3)+k*dh(1,2)+l*dh(1,1) 00995 IF(resp_env%nonperiodic_sys) THEN 00996 not_in_range=.FALSE. 00997 DO m=1,natom 00998 vec_pbc = pbc(r,particles%els(m)%r,cell) 00999 dist(m)=SQRT(SUM(vec_pbc**2)) 01000 CALL get_atomic_kind(atomic_kind=particle_set(m)%atomic_kind,& 01001 kind_number=kind_number) 01002 DO ikind=1,nkind 01003 IF(ikind==kind_number) THEN 01004 rmin=resp_env%rmin_kind(ikind) 01005 rmax=resp_env%rmax_kind(ikind) 01006 EXIT 01007 ENDIF 01008 ENDDO 01009 IF(dist(m)<rmin+1.0E-13_dp) not_in_range(m,1)=.TRUE. 01010 IF(dist(m)>rmax-1.0E-13_dp) not_in_range(m,2)=.TRUE. 01011 ENDDO 01012 IF(ANY(not_in_range(:,1)).OR.ALL(not_in_range(:,2))) CYCLE 01013 ELSE 01014 DO i=1,SIZE(rep_sys) 01015 DO m=1,SIZE(rep_sys(i)%p_resp%atom_surf_list) 01016 in_z=0 01017 in_y=0 01018 in_x=0 01019 iatom=rep_sys(i)%p_resp%atom_surf_list(m) 01020 SELECT CASE(rep_sys(i)%p_resp%my_fit) 01021 CASE(do_resp_x_dir, do_resp_y_dir, do_resp_z_dir) 01022 vec_pbc = pbc(particles%els(iatom)%r,r,cell) 01023 CASE(do_resp_minus_x_dir, do_resp_minus_y_dir, do_resp_minus_z_dir) 01024 vec_pbc = pbc(r,particles%els(iatom)%r,cell) 01025 END SELECT 01026 SELECT CASE(rep_sys(i)%p_resp%my_fit) 01027 !subtract 1.0E-13 to get rid of rounding errors when shifting atoms 01028 CASE(do_resp_x_dir, do_resp_minus_x_dir) 01029 IF(ABS(vec_pbc(3))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_z=1 01030 IF(ABS(vec_pbc(2))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_y=1 01031 IF(vec_pbc(1)>rep_sys(i)%p_resp%range_surf(1)+1.0E-13_dp.AND.& 01032 vec_pbc(1)<rep_sys(i)%p_resp%range_surf(2)-1.0E-13_dp) in_x=1 01033 CASE(do_resp_y_dir,do_resp_minus_y_dir) 01034 IF(ABS(vec_pbc(3))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_z=1 01035 IF(vec_pbc(2)>rep_sys(i)%p_resp%range_surf(1)+1.0E-13_dp.AND.& 01036 vec_pbc(2)<rep_sys(i)%p_resp%range_surf(2)-1.0E-13_dp) in_y=1 01037 IF(ABS(vec_pbc(1))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_x=1 01038 CASE(do_resp_z_dir, do_resp_minus_z_dir) 01039 IF(vec_pbc(3)>rep_sys(i)%p_resp%range_surf(1)+1.0E-13_dp.AND.& 01040 vec_pbc(3)<rep_sys(i)%p_resp%range_surf(2)-1.0E-13_dp) in_z=1 01041 IF(ABS(vec_pbc(2))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_y=1 01042 IF(ABS(vec_pbc(1))<rep_sys(i)%p_resp%length-1.0E-13_dp) in_x=1 01043 END SELECT 01044 IF(in_z*in_y*in_x==1) EXIT 01045 ENDDO 01046 IF(in_z*in_y*in_x==1) EXIT 01047 ENDDO 01048 IF(in_z*in_y*in_x==0) CYCLE 01049 ENDIF 01050 resp_env%npoints=resp_env%npoints+1 01051 IF(resp_env%npoints > now) THEN 01052 now = 2*now 01053 CALL reallocate(resp_env%fitpoints,1,3,1,now) 01054 ENDIF 01055 resp_env%fitpoints(1,resp_env%npoints) = jx 01056 resp_env%fitpoints(2,resp_env%npoints) = jy 01057 resp_env%fitpoints(3,resp_env%npoints) = jz 01058 ENDDO 01059 ENDDO 01060 ENDDO 01061 01062 !print fitting points to file if requested 01063 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01064 resp_section,"PRINT%COORD_FIT_POINTS",error=error),& 01065 cp_p_file))THEN 01066 CALL print_fitting_points(qs_env,resp_env,dh,output_unit,gbo,error) 01067 ENDIF 01068 01069 CALL cp_print_key_finished_output(output_unit,logger,resp_section,& 01070 "PRINT%COORD_FIT_POINTS", error=error) 01071 DEALLOCATE(dist,stat=stat) 01072 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01073 DEALLOCATE(not_in_range,stat=stat) 01074 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01075 01076 CALL timestop(handle) 01077 01078 END SUBROUTINE get_fitting_points 01079 01080 !**************************************************************************** 01081 !\brief calculate vector rhs 01082 !**************************************************************************** 01083 SUBROUTINE calculate_rhs(qs_env,resp_env,rhs,vpot,error) 01084 01085 TYPE(qs_environment_type), POINTER :: qs_env 01086 TYPE(resp_type), POINTER :: resp_env 01087 REAL(KIND=dp), INTENT(INOUT) :: rhs 01088 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vpot 01089 TYPE(cp_error_type), INTENT(inout) :: error 01090 01091 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rhs', 01092 routineP = moduleN//':'//routineN 01093 01094 INTEGER :: handle, ip, jx, jy, jz 01095 LOGICAL :: failure 01096 REAL(KIND=dp) :: dvol, vj 01097 TYPE(pw_type), POINTER :: v_hartree_pw 01098 01099 CALL timeset(routineN,handle) 01100 01101 failure=.FALSE. 01102 01103 NULLIFY(v_hartree_pw) 01104 v_hartree_pw => qs_env%ks_env%v_hartree_rspace%pw 01105 dvol=v_hartree_pw%pw_grid%dvol 01106 !multiply v_hartree and va_rspace and calculate the vector rhs 01107 !taking into account that v_hartree has opposite site; remove v_qmmm 01108 DO ip=1,resp_env%npoints 01109 jx = resp_env%fitpoints(1,ip) 01110 jy = resp_env%fitpoints(2,ip) 01111 jz = resp_env%fitpoints(3,ip) 01112 vj = -v_hartree_pw%cr3d(jx,jy,jz)/dvol 01113 IF (qs_env%qmmm) THEN 01114 !taking into account that v_qmmm has also opposite sign 01115 vj = vj + qs_env%ks_qmmm_env%v_qmmm_rspace%pw%cr3d(jx,jy,jz) 01116 ENDIF 01117 rhs = rhs + 2.0_dp*vj*vpot(ip) 01118 ENDDO 01119 01120 CALL timestop(handle) 01121 01122 END SUBROUTINE calculate_rhs 01123 !**************************************************************************** 01124 !\brief print the atom coordinates and the coordinates of the fitting points 01125 ! to an xyz file 01126 !**************************************************************************** 01127 SUBROUTINE print_fitting_points(qs_env,resp_env,dh,output_unit,gbo,error) 01128 01129 TYPE(qs_environment_type), POINTER :: qs_env 01130 TYPE(resp_type), POINTER :: resp_env 01131 REAL(KIND=dp), INTENT(IN) :: dh(3,3) 01132 INTEGER, INTENT(IN) :: output_unit, gbo(2,3) 01133 TYPE(cp_error_type), INTENT(inout) :: error 01134 01135 CHARACTER(len=*), PARAMETER :: routineN = 'print_fitting_points', 01136 routineP = moduleN//':'//routineN 01137 01138 CHARACTER(LEN=2) :: element_symbol 01139 CHARACTER(LEN=default_path_length) :: filename 01140 INTEGER :: handle, i, iatom, ip, jx, jy, 01141 jz, k, l, my_pos, p, req(6), 01142 stat 01143 INTEGER, DIMENSION(:), POINTER :: tmp_npoints, tmp_size 01144 INTEGER, DIMENSION(:, :), POINTER :: tmp_points 01145 LOGICAL :: failure 01146 REAL(KIND=dp) :: conv, r(3) 01147 TYPE(cp_logger_type), POINTER :: logger 01148 TYPE(cp_para_env_type), POINTER :: para_env 01149 TYPE(particle_type), DIMENSION(:), 01150 POINTER :: particle_set 01151 TYPE(section_vals_type), POINTER :: input, print_key, resp_section 01152 01153 CALL timeset(routineN,handle) 01154 01155 failure=.FALSE. 01156 NULLIFY(para_env,input,logger,resp_section,print_key,particle_set,tmp_size,& 01157 tmp_points,tmp_npoints) 01158 01159 CALL get_qs_env(qs_env, input=input, para_env=para_env,& 01160 particle_set=particle_set,error=error) 01161 01162 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 01163 error=error) 01164 print_key => section_vals_get_subs_vals(resp_section,& 01165 "PRINT%COORD_FIT_POINTS",& 01166 error=error) 01167 logger => cp_error_get_logger(error) 01168 conv=cp_unit_from_cp2k(1.0_dp,"angstrom",error=error) 01169 01170 IF(output_unit>0) THEN 01171 filename=cp_print_key_generate_filename(logger,& 01172 print_key, extension=".xyz",& 01173 my_local=.FALSE.,error=error) 01174 WRITE(unit=output_unit,FMT="(I12,A,/)") SIZE(particle_set),' + nr fit points' 01175 DO iatom=1,SIZE(particle_set) 01176 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01177 element_symbol=element_symbol) 01178 WRITE(UNIT=output_unit,FMT="(A,1X,3F10.5)") element_symbol,& 01179 particle_set(iatom)%r(1:3)*conv 01180 ENDDO 01181 !printing points of proc which is doing the output (should be proc 0) 01182 DO ip=1,resp_env%npoints 01183 jx=resp_env%fitpoints(1,ip) 01184 jy=resp_env%fitpoints(2,ip) 01185 jz=resp_env%fitpoints(3,ip) 01186 l=jx - gbo(1,1) 01187 k=jy - gbo(1,2) 01188 p=jz - gbo(1,3) 01189 r(3)=p*dh(3,3)+k*dh(3,2)+l*dh(3,1) 01190 r(2)=p*dh(2,3)+k*dh(2,2)+l*dh(2,1) 01191 r(1)=p*dh(1,3)+k*dh(1,2)+l*dh(1,1) 01192 r(:)=r(:)*conv 01193 WRITE(UNIT=output_unit,FMT="(A,2X,3F10.5)") "X", r(1), r(2), r(3) 01194 ENDDO 01195 ENDIF 01196 01197 ALLOCATE(tmp_size(1),stat=stat) 01198 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 01199 ALLOCATE(tmp_npoints(1),stat=stat) 01200 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 01201 01202 !sending data of all other prosc to proc which makes the output (proc 0) 01203 IF(output_unit>0) THEN 01204 my_pos=para_env%mepos 01205 DO i=1,para_env%num_pe 01206 IF(my_pos==i-1) CYCLE 01207 CALL mp_irecv(msgout=tmp_size,source=i-1,comm=para_env%group,& 01208 request=req(1)) 01209 CALL mp_wait(req(1)) 01210 ALLOCATE(tmp_points(3,tmp_size(1)),stat=stat) 01211 CPPostcondition(stat==0,cp_failure_level,routineP,error,Failure) 01212 CALL mp_irecv(msgout=tmp_points,source=i-1,comm=para_env%group,& 01213 request=req(3)) 01214 CALL mp_wait(req(3)) 01215 CALL mp_irecv(msgout=tmp_npoints,source=i-1,comm=para_env%group,& 01216 request=req(5)) 01217 CALL mp_wait(req(5)) 01218 DO ip=1,tmp_npoints(1) 01219 jx=tmp_points(1,ip) 01220 jy=tmp_points(2,ip) 01221 jz=tmp_points(3,ip) 01222 l=jx - gbo(1,1) 01223 k=jy - gbo(1,2) 01224 p=jz - gbo(1,3) 01225 r(3)=p*dh(3,3)+k*dh(3,2)+l*dh(3,1) 01226 r(2)=p*dh(2,3)+k*dh(2,2)+l*dh(2,1) 01227 r(1)=p*dh(1,3)+k*dh(1,2)+l*dh(1,1) 01228 r(:)=r(:)*conv 01229 WRITE(UNIT=output_unit,FMT="(A,2X,3F10.5)") "X", r(1), r(2),r(3) 01230 ENDDO 01231 DEALLOCATE(tmp_points,stat=stat) 01232 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01233 ENDDO 01234 ELSE 01235 tmp_size(1)=SIZE(resp_env%fitpoints,2) 01236 !para_env%source should be 0 01237 CALL mp_isend(msgin=tmp_size,dest=para_env%source,comm=para_env%group,& 01238 request=req(2)) 01239 CALL mp_wait(req(2)) 01240 CALL mp_isend(msgin=resp_env%fitpoints,dest=para_env%source,comm=para_env%group,& 01241 request=req(4)) 01242 CALL mp_wait(req(4)) 01243 tmp_npoints(1)=resp_env%npoints 01244 CALL mp_isend(msgin=tmp_npoints,dest=para_env%source,comm=para_env%group,& 01245 request=req(6)) 01246 CALL mp_wait(req(6)) 01247 ENDIF 01248 01249 DEALLOCATE(tmp_size,stat=stat) 01250 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01251 DEALLOCATE(tmp_npoints,stat=stat) 01252 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01253 01254 CALL timestop(handle) 01255 01256 END SUBROUTINE print_fitting_points 01257 !**************************************************************************** 01258 !\brief add restraints and constraints 01259 !**************************************************************************** 01260 SUBROUTINE add_restraints_and_constraints(qs_env,resp_env,rest_section,& 01261 subsys,natom,cons_section,particle_set,error) 01262 01263 TYPE(qs_environment_type), POINTER :: qs_env 01264 TYPE(resp_type), POINTER :: resp_env 01265 TYPE(section_vals_type), POINTER :: rest_section 01266 TYPE(cp_subsys_type), POINTER :: subsys 01267 INTEGER, INTENT(IN) :: natom 01268 TYPE(section_vals_type), POINTER :: cons_section 01269 TYPE(particle_type), DIMENSION(:), 01270 POINTER :: particle_set 01271 TYPE(cp_error_type), INTENT(inout) :: error 01272 01273 CHARACTER(len=*), PARAMETER :: 01274 routineN = 'add_restraints_and_constraints', 01275 routineP = moduleN//':'//routineN 01276 01277 INTEGER :: handle, i, k, m, ncons_v, 01278 stat, z 01279 INTEGER, DIMENSION(:), POINTER :: atom_list_cons, atom_list_res 01280 LOGICAL :: explicit_coeff, failure 01281 REAL(KIND=dp) :: my_atom_coef(2), strength, 01282 TARGET 01283 REAL(KIND=dp), DIMENSION(:), POINTER :: atom_coef 01284 01285 CALL timeset(routineN,handle) 01286 01287 NULLIFY(atom_coef,atom_list_res,atom_list_cons) 01288 failure=.FALSE. 01289 01290 ! add the restraints 01291 DO i=1,resp_env%nrest_sec 01292 CALL section_vals_val_get(rest_section,"TARGET",i_rep_section=i,r_val=TARGET,error=error) 01293 CALL section_vals_val_get(rest_section,"STRENGTH",i_rep_section=i,r_val=strength,error=error) 01294 CALL build_atom_list(rest_section,subsys,atom_list_res,i,error) 01295 CALL section_vals_val_get(rest_section,"ATOM_COEF",i_rep_section=i,explicit=explicit_coeff,& 01296 error=error) 01297 IF (explicit_coeff) THEN 01298 CALL section_vals_val_get(rest_section,"ATOM_COEF",i_rep_section=i,r_vals=atom_coef,& 01299 error=error) 01300 CPPrecondition(SIZE(atom_list_res)==SIZE(atom_coef),cp_failure_level,routineP,error,failure) 01301 ENDIF 01302 DO m=1,SIZE(atom_list_res) 01303 IF (explicit_coeff) THEN 01304 DO k=1,SIZE(atom_list_res) 01305 resp_env%matrix(atom_list_res(m),atom_list_res(k))=& 01306 resp_env%matrix(atom_list_res(m),atom_list_res(k))+ & 01307 atom_coef(m)*atom_coef(k)*2.0_dp*strength 01308 ENDDO 01309 resp_env%rhs(atom_list_res(m))=resp_env%rhs(atom_list_res(m))+& 01310 2.0_dp*TARGET*strength*atom_coef(m) 01311 ELSE 01312 resp_env%matrix(atom_list_res(m),atom_list_res(m))=& 01313 resp_env%matrix(atom_list_res(m),atom_list_res(m))+ & 01314 2.0_dp*strength 01315 resp_env%rhs(atom_list_res(m))=resp_env%rhs(atom_list_res(m))+& 01316 2.0_dp*TARGET*strength 01317 ENDIF 01318 ENDDO 01319 DEALLOCATE(atom_list_res,stat=stat) 01320 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01321 ENDDO 01322 01323 ! if heavies are restrained to zero, add these as well 01324 IF (resp_env%rheavies) THEN 01325 DO i=1,natom 01326 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind,z=z) 01327 IF (z.NE.1) THEN 01328 resp_env%matrix(i,i)=resp_env%matrix(i,i)+2.0_dp*resp_env%rheavies_strength 01329 ENDIF 01330 ENDDO 01331 ENDIF 01332 01333 ! add the constraints 01334 ncons_v=0 01335 ncons_v=ncons_v+natom 01336 IF (resp_env%itc) THEN 01337 ncons_v=ncons_v+1 01338 resp_env%matrix(1:natom,ncons_v)=1.0_dp 01339 resp_env%matrix(ncons_v,1:natom)=1.0_dp 01340 resp_env%rhs(ncons_v)=qs_env%dft_control%charge 01341 ENDIF 01342 01343 DO i=1,resp_env%ncons_sec 01344 CALL build_atom_list(cons_section,subsys,atom_list_cons,i,error) 01345 IF(.NOT.resp_env%equal_charges) THEN 01346 ncons_v=ncons_v+1 01347 CALL section_vals_val_get(cons_section,"ATOM_COEF",i_rep_section=i,r_vals=atom_coef,error=error) 01348 CALL section_vals_val_get(cons_section,"TARGET",i_rep_section=i,r_val=TARGET,error=error) 01349 CPPrecondition(SIZE(atom_list_cons)==SIZE(atom_coef),cp_failure_level,routineP,error,failure) 01350 DO m=1,SIZE(atom_list_cons) 01351 resp_env%matrix(atom_list_cons(m),ncons_v)=atom_coef(m) 01352 resp_env%matrix(ncons_v,atom_list_cons(m))=atom_coef(m) 01353 ENDDO 01354 resp_env%rhs(ncons_v)=TARGET 01355 ELSE 01356 my_atom_coef(1)=1.0_dp 01357 my_atom_coef(2)=-1.0_dp 01358 DO k=2,SIZE(atom_list_cons) 01359 ncons_v=ncons_v+1 01360 resp_env%matrix(atom_list_cons(1),ncons_v)=my_atom_coef(1) 01361 resp_env%matrix(ncons_v,atom_list_cons(1))=my_atom_coef(1) 01362 resp_env%matrix(atom_list_cons(k),ncons_v)=my_atom_coef(2) 01363 resp_env%matrix(ncons_v,atom_list_cons(k))=my_atom_coef(2) 01364 resp_env%rhs(ncons_v)=0.0_dp 01365 ENDDO 01366 ENDIF 01367 DEALLOCATE(atom_list_cons,stat=stat) 01368 CPPostconditionNoFail(stat==0,cp_failure_level,routineP,error) 01369 ENDDO 01370 CALL timestop(handle) 01371 01372 END SUBROUTINE add_restraints_and_constraints 01373 !**************************************************************************** 01374 !\brief print input information 01375 !**************************************************************************** 01376 SUBROUTINE print_resp_parameter_info(qs_env,resp_env,rep_sys,my_per,error) 01377 01378 TYPE(qs_environment_type), POINTER :: qs_env 01379 TYPE(resp_type), POINTER :: resp_env 01380 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys 01381 INTEGER, INTENT(IN) :: my_per 01382 TYPE(cp_error_type), INTENT(inout) :: error 01383 01384 CHARACTER(len=*), PARAMETER :: routineN = 'print_resp_parameter_info', 01385 routineP = moduleN//':'//routineN 01386 01387 INTEGER :: handle, i, output_unit 01388 REAL(KIND=dp) :: conv, eta_conv 01389 TYPE(cp_logger_type), POINTER :: logger 01390 TYPE(section_vals_type), POINTER :: input, resp_section 01391 01392 CALL timeset(routineN,handle) 01393 NULLIFY(logger,input,resp_section) 01394 01395 CALL get_qs_env(qs_env, input=input, error=error) 01396 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 01397 error=error) 01398 logger => cp_error_get_logger(error) 01399 output_unit=cp_print_key_unit_nr(logger,resp_section,"PRINT%PROGRAM_RUN_INFO",& 01400 extension=".resp",error=error) 01401 01402 conv=cp_unit_from_cp2k(1.0_dp,"angstrom",error=error) 01403 IF(.NOT.my_per==use_perd_none) THEN 01404 eta_conv=cp_unit_from_cp2k(resp_env%eta,"angstrom",power=-2,error=error) 01405 ENDIF 01406 01407 IF (output_unit>0) THEN 01408 WRITE(output_unit,'(/,1X,A,/)') "STARTING RESP FIT" 01409 IF(.NOT.resp_env%equal_charges) THEN 01410 WRITE(output_unit,'(T3,A,T75,I6)') "Number of explicit constraints: ",resp_env%ncons_sec 01411 ELSE 01412 IF(resp_env%itc) THEN 01413 WRITE(output_unit,'(T3,A,T75,I6)') "Number of explicit constraints: ",resp_env%ncons-1 01414 ELSE 01415 WRITE(output_unit,'(T3,A,T75,I6)') "Number of explicit constraints: ",resp_env%ncons 01416 ENDIF 01417 ENDIF 01418 WRITE(output_unit,'(T3,A,T75,I6)') "Number of explicit restraints: ",resp_env%nrest_sec 01419 WRITE(output_unit,'(T3,A,T80,A)') "Constrain total charge ",MERGE("T","F",resp_env%itc) 01420 WRITE(output_unit,'(T3,A,T80,A)') "Restrain heavy atoms ",MERGE("T","F",resp_env%rheavies) 01421 IF (resp_env%rheavies) THEN 01422 WRITE(output_unit,'(T3,A,T71,F10.6)') "Heavy atom restraint strength: ",& 01423 resp_env%rheavies_strength 01424 ENDIF 01425 WRITE(output_unit,'(T3,A,T66,3I5)') "Stride: ",resp_env%stride 01426 IF(resp_env%nonperiodic_sys) THEN 01427 IF(ALL(resp_env%rmin_kind==resp_env%rmin_kind(1)))& 01428 WRITE(output_unit,'(T3,A,T71,F10.5)') "Rmin [angstrom]: ",resp_env%rmin_kind(1)*conv 01429 IF(ALL(resp_env%rmax_kind==resp_env%rmax_kind(1)))& 01430 WRITE(output_unit,'(T3,A,T71,F10.5)') "Rmax [angstrom]: ",resp_env%rmax_kind(1)*conv 01431 WRITE(output_unit,'(T3,A,T51,3F10.5)') "Box min [angstrom]: ",resp_env%box_low(1:3)*conv 01432 WRITE(output_unit,'(T3,A,T51,3F10.5)') "Box max [angstrom]: ",resp_env%box_hi(1:3)*conv 01433 ELSE 01434 WRITE(output_unit,'(2X,A,F10.5)') "Index of atoms defining the surface: " 01435 DO i=1, SIZE(rep_sys) 01436 IF(i>1.AND.ALL(rep_sys(i)%p_resp%atom_surf_list==rep_sys(1)%p_resp%atom_surf_list)) EXIT 01437 WRITE(output_unit,'(7X,10I6)') rep_sys(i)%p_resp%atom_surf_list 01438 ENDDO 01439 DO i=1, SIZE(rep_sys) 01440 IF(i>1.AND.ALL(rep_sys(i)%p_resp%range_surf==rep_sys(1)%p_resp%range_surf)) EXIT 01441 WRITE(output_unit,'(T3,A,T61,2F10.5)')& 01442 "Range for sampling above the surface [angstrom]:",& 01443 rep_sys(i)%p_resp%range_surf(1:2)*conv 01444 ENDDO 01445 DO i=1, SIZE(rep_sys) 01446 IF(i>1.AND.rep_sys(i)%p_resp%length==rep_sys(1)%p_resp%length) EXIT 01447 WRITE(output_unit,'(T3,A,T71,F10.5)') "Length of sampling box above each"//& 01448 " surface atom [angstrom]: ",rep_sys(i)%p_resp%length*conv 01449 ENDDO 01450 ENDIF 01451 IF(.NOT.my_per==use_perd_none) THEN 01452 WRITE(output_unit,'(T3,A,T71,F10.5)') "Width of Gaussian charge"//& 01453 " distribution [angstrom^-2]: ", eta_conv 01454 ENDIF 01455 CALL m_flush(output_unit) 01456 ENDIF 01457 CALL cp_print_key_finished_output(output_unit,logger,resp_section,& 01458 "PRINT%PROGRAM_RUN_INFO", error=error) 01459 01460 CALL timestop(handle) 01461 01462 END SUBROUTINE print_resp_parameter_info 01463 !**************************************************************************** 01464 !\brief print RESP charges to an extra file or to the normal output file 01465 !**************************************************************************** 01466 SUBROUTINE print_resp_charges(qs_env,resp_env,output_runinfo,natom,error) 01467 01468 TYPE(qs_environment_type), POINTER :: qs_env 01469 TYPE(resp_type), POINTER :: resp_env 01470 INTEGER, INTENT(IN) :: output_runinfo, natom 01471 TYPE(cp_error_type), INTENT(inout) :: error 01472 01473 CHARACTER(len=*), PARAMETER :: routineN = 'print_resp_charges', 01474 routineP = moduleN//':'//routineN 01475 01476 CHARACTER(LEN=default_path_length) :: filename 01477 INTEGER :: handle, output_file 01478 LOGICAL :: failure 01479 TYPE(cp_logger_type), POINTER :: logger 01480 TYPE(particle_type), DIMENSION(:), 01481 POINTER :: particle_set 01482 TYPE(section_vals_type), POINTER :: input, print_key, resp_section 01483 01484 CALL timeset(routineN,handle) 01485 01486 failure=.FALSE. 01487 NULLIFY(particle_set,input,logger,resp_section,print_key) 01488 01489 CALL get_qs_env(qs_env,input=input,particle_set=particle_set,error=error) 01490 01491 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 01492 error=error) 01493 print_key => section_vals_get_subs_vals(resp_section,& 01494 "PRINT%RESP_CHARGES_TO_FILE",& 01495 error=error) 01496 logger => cp_error_get_logger(error) 01497 01498 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01499 resp_section,"PRINT%RESP_CHARGES_TO_FILE",error=error),& 01500 cp_p_file)) THEN 01501 output_file=cp_print_key_unit_nr(logger,resp_section,& 01502 "PRINT%RESP_CHARGES_TO_FILE",& 01503 extension=".resp",& 01504 file_status="REPLACE",& 01505 file_action="WRITE",& 01506 file_form="FORMATTED",& 01507 error=error) 01508 IF(output_file>0) THEN 01509 filename = cp_print_key_generate_filename(logger,& 01510 print_key, extension=".resp", & 01511 my_local=.FALSE.,error=error) 01512 CALL print_atomic_charges(particle_set,output_file,title="RESP charges:",& 01513 atomic_charges=resp_env%rhs(1:natom)) 01514 IF(output_runinfo>0) WRITE(output_runinfo,'(2X,A,/)') "PRINTED RESP CHARGES TO FILE" 01515 ENDIF 01516 01517 CALL cp_print_key_finished_output(output_file,logger,resp_section,& 01518 "PRINT%RESP_CHARGES_TO_FILE", error=error) 01519 ELSE 01520 CALL print_atomic_charges(particle_set,output_runinfo,title="RESP charges:",& 01521 atomic_charges=resp_env%rhs(1:natom)) 01522 ENDIF 01523 01524 CALL timestop(handle) 01525 END SUBROUTINE print_resp_charges 01526 !**************************************************************************** 01527 !\brief print potential generates by RESP charges to file 01528 !**************************************************************************** 01529 SUBROUTINE print_pot_from_resp_charges(qs_env,resp_env,particles,natom,output_runinfo,error) 01530 01531 TYPE(qs_environment_type), POINTER :: qs_env 01532 TYPE(resp_type), POINTER :: resp_env 01533 TYPE(particle_list_type), POINTER :: particles 01534 INTEGER, INTENT(IN) :: natom, output_runinfo 01535 TYPE(cp_error_type), INTENT(inout) :: error 01536 01537 CHARACTER(len=*), PARAMETER :: routineN = 'print_pot_from_resp_charges', 01538 routineP = moduleN//':'//routineN 01539 01540 CHARACTER(LEN=default_path_length) :: my_pos_cube 01541 INTEGER :: handle, ip, jx, jy, jz, 01542 unit_nr 01543 LOGICAL :: append_cube, failure 01544 REAL(KIND=dp) :: normalize_factor, rms, rrms, 01545 sum_diff, sum_hartree, udvol 01546 TYPE(cp_logger_type), POINTER :: logger 01547 TYPE(cp_para_env_type), POINTER :: para_env 01548 TYPE(pw_env_type), POINTER :: pw_env 01549 TYPE(pw_p_type) :: aux_r, rho_resp, 01550 v_resp_gspace, v_resp_rspace 01551 TYPE(pw_poisson_type), POINTER :: poisson_env 01552 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01553 TYPE(section_vals_type), POINTER :: input, print_key, resp_section 01554 01555 CALL timeset(routineN,handle) 01556 01557 failure=.FALSE. 01558 NULLIFY(auxbas_pw_pool,logger,pw_env,poisson_env,input,print_key,& 01559 para_env,resp_section) 01560 normalize_factor=SQRT((resp_env%eta/pi)**3) 01561 CALL get_qs_env(qs_env=qs_env,input=input,para_env=para_env,& 01562 pw_env=pw_env,error=error) 01563 resp_section => section_vals_get_subs_vals(input,"PROPERTIES%RESP",& 01564 error=error) 01565 print_key => section_vals_get_subs_vals(resp_section,& 01566 "PRINT%V_RESP_CUBE",error=error) 01567 logger => cp_error_get_logger(error) 01568 IF (BTEST(cp_print_key_should_output(logger%iter_info,resp_section,& 01569 "PRINT%V_RESP_CUBE",error=error),cp_p_file)) THEN 01570 ! calculate potential generated from RESP charges 01571 CALL pw_env_get(pw_env,auxbas_pw_pool=auxbas_pw_pool,& 01572 poisson_env=poisson_env,error=error) 01573 01574 CALL pw_pool_create_pw(auxbas_pw_pool,& 01575 rho_resp%pw,& 01576 use_data=COMPLEXDATA1D,& 01577 in_space=RECIPROCALSPACE,& 01578 error=error) 01579 CALL pw_pool_create_pw(auxbas_pw_pool,& 01580 v_resp_gspace%pw, & 01581 use_data=COMPLEXDATA1D,& 01582 in_space=RECIPROCALSPACE,& 01583 error=error) 01584 CALL pw_pool_create_pw(auxbas_pw_pool,& 01585 v_resp_rspace%pw,& 01586 use_data=REALDATA3D,& 01587 in_space=REALSPACE,& 01588 error=error) 01589 01590 CALL pw_zero(rho_resp%pw,error=error) 01591 CALL calculate_rho_resp_all(rho_resp,resp_env%rhs,natom,& 01592 resp_env%eta,qs_env,error) 01593 CALL pw_zero(v_resp_gspace%pw, error=error) 01594 CALL pw_poisson_solve(poisson_env,rho_resp%pw,& 01595 vhartree=v_resp_gspace%pw,error=error) 01596 CALL pw_zero(v_resp_rspace%pw, error=error) 01597 CALL pw_transfer(v_resp_gspace%pw,v_resp_rspace%pw,error=error) 01598 CALL pw_scale(v_resp_rspace%pw,v_resp_rspace%pw%pw_grid%dvol,& 01599 error=error) 01600 CALL pw_scale(v_resp_rspace%pw,-normalize_factor,& 01601 error=error) 01602 CALL pw_release(v_resp_gspace%pw,error=error) 01603 CALL pw_release(rho_resp%pw,error=error) 01604 !now print the v_resp_rspace%pw to a cube file 01605 CALL pw_pool_create_pw(auxbas_pw_pool,aux_r%pw,& 01606 use_data = REALDATA3D,& 01607 in_space = REALSPACE, error=error) 01608 append_cube=section_get_lval(resp_section,"PRINT%V_RESP_CUBE%APPEND",error=error) 01609 my_pos_cube="REWIND" 01610 IF(append_cube) THEN 01611 my_pos_cube="APPEND" 01612 END IF 01613 unit_nr=cp_print_key_unit_nr(logger,resp_section,& 01614 "PRINT%V_RESP_CUBE",& 01615 extension=".cube",& 01616 file_position=my_pos_cube,error=error) 01617 udvol = 1.0_dp/v_resp_rspace%pw%pw_grid%dvol 01618 CALL pw_copy(v_resp_rspace%pw,aux_r%pw, error=error) 01619 CALL pw_scale(aux_r%pw,udvol,error=error) 01620 CALL pw_to_cube(aux_r%pw,unit_nr,"RESP POTENTIAL",particles=particles,& 01621 stride=section_get_ivals(resp_section,& 01622 "PRINT%V_RESP_CUBE%STRIDE",error=error),& 01623 error=error) 01624 CALL cp_print_key_finished_output(unit_nr,logger,resp_section,& 01625 "PRINT%V_RESP_CUBE",error=error) 01626 CALL pw_pool_give_back_pw(auxbas_pw_pool,aux_r%pw, error=error) 01627 01628 !RMS and RRMS 01629 sum_diff=0.0_dp 01630 sum_hartree=0.0_dp 01631 rms=0.0_dp 01632 rrms=0.0_dp 01633 DO ip=1,resp_env%npoints_proc 01634 jx = resp_env%fitpoints(1,ip) 01635 jy = resp_env%fitpoints(2,ip) 01636 jz = resp_env%fitpoints(3,ip) 01637 sum_diff=sum_diff+(qs_env%ks_env%v_hartree_rspace%pw%cr3d(jx,jy,jz)-& 01638 v_resp_rspace%pw%cr3d(jx,jy,jz))**2 01639 sum_hartree=sum_hartree+qs_env%ks_env%v_hartree_rspace%pw%cr3d(jx,jy,jz)**2 01640 ENDDO 01641 CALL mp_sum(sum_diff,para_env%group) 01642 CALL mp_sum(sum_hartree,para_env%group) 01643 rms=SQRT(sum_diff/resp_env%npoints) 01644 rrms=SQRT(sum_diff/sum_hartree) 01645 IF(output_runinfo > 0) THEN 01646 WRITE(output_runinfo,'(2X,A,T69,ES12.5)') "Root-mean-square (RMS) "//& 01647 "error of RESP fit:",rms 01648 WRITE(output_runinfo,'(2X,A,T69,ES12.5,/)') "Relative root-mean-square "//& 01649 "(RRMS) error of RESP fit:",rrms 01650 ENDIF 01651 CALL pw_release(v_resp_rspace%pw,error=error) 01652 ENDIF 01653 01654 CALL timestop(handle) 01655 END SUBROUTINE print_pot_from_resp_charges 01656 END MODULE qs_resp
1.7.3