|
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 ! ***************************************************************************** 00015 MODULE force_env_methods 00016 00017 USE atprop_types, ONLY: atprop_create,& 00018 atprop_init,& 00019 atprop_type 00020 USE cell_types, ONLY: cell_clone,& 00021 cell_create,& 00022 cell_release,& 00023 cell_type,& 00024 compare_cells,& 00025 init_cell,& 00026 real_to_scaled,& 00027 scaled_to_real 00028 USE constraint_fxd, ONLY: fix_atom_control 00029 USE constraint_vsite, ONLY: vsite_force_control 00030 USE cp_iter_types, ONLY: cp_iteration_info_copy_iter 00031 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00032 cp_print_key_unit_nr 00033 USE cp_para_env, ONLY: cp_para_env_retain 00034 USE cp_para_types, ONLY: cp_para_env_type 00035 USE cp_result_methods, ONLY: cp_results_erase,& 00036 cp_results_mp_bcast,& 00037 get_results 00038 USE cp_result_types, ONLY: cp_result_copy,& 00039 cp_result_create,& 00040 cp_result_p_type,& 00041 cp_result_release,& 00042 cp_result_type 00043 USE cp_subsys_types, ONLY: cp_subsys_get,& 00044 cp_subsys_p_type,& 00045 cp_subsys_retain,& 00046 cp_subsys_type 00047 USE cpot_types, ONLY: cpot_calc 00048 USE eip_environment_types, ONLY: eip_env_get,& 00049 eip_env_retain,& 00050 eip_environment_type 00051 USE eip_silicon, ONLY: eip_bazant,& 00052 eip_lenosky 00053 USE ep_types, ONLY: ep_env_calc_e_f,& 00054 ep_env_create,& 00055 ep_env_release,& 00056 ep_env_retain,& 00057 ep_env_type 00058 USE external_potential_methods, ONLY: add_external_potential 00059 USE f77_blas 00060 USE fist_environment_types, ONLY: fist_env_get,& 00061 fist_env_retain,& 00062 fist_environment_type 00063 USE fist_force, ONLY: fist_force_control 00064 USE force_env_types, ONLY: & 00065 force_env_get, force_env_get_natom, force_env_get_pos, & 00066 force_env_p_type, force_env_set, force_env_set_cell, force_env_type, & 00067 use_eip_force, use_ep_force, use_fist_force, use_mixed_force, & 00068 use_prog_name, use_qmmm, use_qs_force 00069 USE force_env_utils, ONLY: rescale_forces,& 00070 write_forces,& 00071 write_stress_tensor 00072 USE force_fields_util, ONLY: get_generic_info 00073 USE fp_methods, ONLY: fp_eval 00074 USE fparser, ONLY: EvalErrType,& 00075 evalf,& 00076 evalfd,& 00077 finalizef,& 00078 initf,& 00079 parsef 00080 USE global_types, ONLY: global_environment_type,& 00081 globenv_retain 00082 USE input_constants, ONLY: & 00083 debug_run, do_fm_mom_conserv_QM, do_fm_mom_conserv_buffer, & 00084 do_fm_mom_conserv_core, do_fm_mom_conserv_equal_a, & 00085 do_fm_mom_conserv_equal_f, do_fm_mom_conserv_none, & 00086 do_stress_analytical, do_stress_diagonal_anal, & 00087 do_stress_diagonal_numer, do_stress_none, do_stress_numerical, & 00088 low_print_level, mix_coupled, mix_generic, mix_linear_combination, & 00089 mix_minimum, mix_restrained, use_bazant_eip, use_lenosky_eip 00090 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00091 section_vals_retain,& 00092 section_vals_type,& 00093 section_vals_val_get,& 00094 section_vals_val_set 00095 USE kinds, ONLY: default_path_length,& 00096 default_string_length,& 00097 dp 00098 USE message_passing, ONLY: mp_sum,& 00099 mp_sync 00100 USE metadynamics_types, ONLY: meta_env_retain,& 00101 meta_env_type 00102 USE mixed_energy_types, ONLY: mixed_energy_type,& 00103 mixed_force_type 00104 USE mixed_environment_types, ONLY: get_mixed_env,& 00105 mixed_env_retain,& 00106 mixed_environment_type 00107 USE mixed_environment_utils, ONLY: get_subsys_map_index,& 00108 mixed_map_forces 00109 USE particle_list_types, ONLY: particle_list_p_type,& 00110 particle_list_type 00111 USE particle_types, ONLY: particle_type 00112 USE physcon, ONLY: debye 00113 USE qmmm_gpw_energy, ONLY: qmmm_el_coupling 00114 USE qmmm_gpw_forces, ONLY: qmmm_forces 00115 USE qmmm_links_methods, ONLY: qmmm_added_chrg_coord,& 00116 qmmm_added_chrg_forces,& 00117 qmmm_link_Imomm_coord,& 00118 qmmm_link_Imomm_forces 00119 USE qmmm_types, ONLY: & 00120 fist_subsys, force_mixing_core_subsys, force_mixing_extended_subsys, & 00121 force_mixing_label_QM_core, force_mixing_label_QM_dynamics, & 00122 force_mixing_label_buffer, primary_subsys, qmmm_env_qm_retain, & 00123 qmmm_env_qm_type, qmmm_links_type, qs_subsys 00124 USE qmmm_util, ONLY: apply_qmmm_translate,& 00125 apply_qmmm_walls,& 00126 qmmm_force_mixing_active 00127 USE qs_energy, ONLY: qs_energies 00128 USE qs_environment_types, ONLY: get_qs_env,& 00129 qs_env_retain,& 00130 qs_environment_type,& 00131 set_qs_env 00132 USE qs_force, ONLY: qs_forces 00133 USE qs_ks_qmmm_methods, ONLY: ks_qmmm_env_rebuild 00134 USE restraint, ONLY: restraint_control 00135 USE se_ga_tools, ONLY: se_ga_release 00136 USE string_utilities, ONLY: compress 00137 USE virial_types, ONLY: & 00138 cp_virial, sym_virial, virial_create, virial_p_type, virial_release, & 00139 virial_retain, virial_set, virial_type, zero_virial 00140 #include "cp_common_uses.h" 00141 00142 IMPLICIT NONE 00143 00144 PRIVATE 00145 00146 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_methods' 00147 00148 PUBLIC :: force_env_create,& 00149 ep_create_force_env,& 00150 force_env_calc_energy_force,& 00151 force_env_calc_num_pressure 00152 00153 INTEGER, SAVE, PRIVATE :: last_force_env_id=0 00154 00155 CONTAINS 00156 00157 ! ***************************************************************************** 00167 RECURSIVE SUBROUTINE force_env_calc_energy_force ( force_env, calc_force, & 00168 consistent_energies, skip_external_control, eval_energy_forces, & 00169 require_consistent_energy_force, error) 00170 00171 TYPE(force_env_type), POINTER :: force_env 00172 LOGICAL, INTENT(IN), OPTIONAL :: calc_force, consistent_energies, 00173 skip_external_control, eval_energy_forces, 00174 require_consistent_energy_force 00175 TYPE(cp_error_type), INTENT(inout) :: error 00176 00177 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_calc_energy_force', 00178 routineP = moduleN//':'//routineN 00179 REAL(kind=dp), PARAMETER :: ateps = 1.0E-6_dp 00180 00181 INTEGER :: i, j, nat, ndigits, 00182 output_unit, print_forces, 00183 stat 00184 LOGICAL :: calculate_forces, 00185 energy_consistency, eval_ef, 00186 failure, my_skip 00187 REAL(KIND=dp) :: checksum, e_pot, sum_energy, 00188 sum_pv_virial, 00189 sum_stress_tensor 00190 REAL(KIND=dp), DIMENSION(3) :: grand_total_force, total_force 00191 REAL(KIND=dp), DIMENSION(3, 3) :: atomic_stress_tensor, 00192 diff_stress_tensor 00193 REAL(KIND=dp), DIMENSION(:), POINTER :: pos 00194 TYPE(atprop_type), POINTER :: atprop_env 00195 TYPE(cell_type), POINTER :: cell 00196 TYPE(cp_logger_type), POINTER :: logger 00197 TYPE(cp_result_type), POINTER :: results 00198 TYPE(cp_subsys_type), POINTER :: subsys 00199 TYPE(particle_list_type), POINTER :: core_particles, particles, 00200 shell_particles 00201 TYPE(virial_type), POINTER :: virial 00202 00203 NULLIFY (logger,results) 00204 logger => cp_error_get_logger(error) 00205 failure = .FALSE. 00206 eval_ef = .TRUE. 00207 my_skip = .FALSE. 00208 calculate_forces = .TRUE. 00209 energy_consistency = .FALSE. 00210 IF (PRESENT(eval_energy_forces)) eval_ef = eval_energy_forces 00211 IF (PRESENT(skip_external_control)) my_skip = skip_external_control 00212 IF (PRESENT(calc_force)) calculate_forces = calc_force 00213 IF (PRESENT(consistent_energies)) energy_consistency = consistent_energies 00214 00215 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00216 IF (.NOT.failure) THEN 00217 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00218 END IF 00219 00220 CALL force_env_get(force_env,virial=virial,error=error) 00221 CALL force_env_set(force_env,additional_potential=0.0_dp,error=error) 00222 IF (virial%pv_availability) CALL zero_virial(virial,reset=.FALSE.) 00223 00224 NULLIFY (atprop_env) 00225 CALL force_env_get(force_env,atprop_env=atprop_env,error=error) 00226 nat=force_env_get_natom(force_env,error=error) 00227 CALL atprop_init(atprop_env,nat,error) 00228 00229 IF (.NOT.failure) THEN 00230 IF (eval_ef) THEN 00231 SELECT CASE ( force_env%in_use ) 00232 CASE ( use_fist_force ) 00233 CALL fist_force_control( force_env%fist_env, virial, atprop_env, force_env%para_env, & 00234 force_env_section=force_env%force_env_section, error=error) 00235 CALL fist_env_get(fist_env=force_env%fist_env,results=results,error=error) 00236 CASE (use_ep_force) 00237 CALL ep_env_calc_e_f(force_env%ep_env,calculate_forces,error=error) 00238 CASE ( use_qs_force ) 00239 CALL set_qs_env(qs_env=force_env%qs_env,atprop=atprop_env,error=error) 00240 IF (.NOT.calculate_forces) THEN 00241 CALL qs_energies(qs_env=force_env%qs_env, consistent_energies=energy_consistency, & 00242 calc_forces=calculate_forces, error=error) 00243 ELSE 00244 CALL qs_forces(force_env%qs_env, force_env%globenv, error=error) 00245 END IF 00246 ! GA option 00247 CALL se_ga_release ( force_env%qs_env, error ) 00248 CALL get_qs_env(qs_env=force_env%qs_env,results=results,error=error) 00249 CASE (use_eip_force) 00250 IF (force_env%eip_env%eip_model == use_lenosky_eip) THEN 00251 CALL eip_lenosky(force_env, error=error) 00252 ELSE IF (force_env%eip_env%eip_model == use_bazant_eip) THEN 00253 CALL eip_bazant(force_env, error=error) 00254 END IF 00255 CASE ( use_qmmm ) 00256 CALL qmmm_energy_and_forces(force_env,calculate_forces,require_consistent_energy_force,& 00257 error=error) 00258 CASE ( use_mixed_force ) 00259 CALL mixed_energy_forces(force_env,calculate_forces,error=error) 00260 CASE default 00261 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00262 END SELECT 00263 END IF 00264 ! In case it is requested, we evaluate the stress tensor numerically 00265 IF (virial%pv_availability.AND.calculate_forces) THEN 00266 IF (virial%pv_numer) THEN 00267 ! Compute the numerical stress tensor 00268 CALL force_env_calc_num_pressure(force_env, error=error) 00269 ELSE 00270 ! Symmetrize analytical stress tensor 00271 CALL sym_virial(virial, error) 00272 END IF 00273 END IF 00274 ! Some additional tasks.. 00275 IF (.NOT.my_skip) THEN 00276 ! Flexible Partitioning 00277 IF (ASSOCIATED(force_env%fp_env)) THEN 00278 IF (force_env%fp_env%use_fp) THEN 00279 CALL force_env_get(force_env,cell=cell,error=error) 00280 CALL fp_eval(force_env%fp_env,force_env%subsys,cell,error=error) 00281 ENDIF 00282 ENDIF 00283 ! Constraints ONLY of Fixed Atom type 00284 CALL fix_atom_control(force_env, error=error) 00285 ! All Restraints 00286 CALL restraint_control(force_env, error=error) 00287 ! Virtual Sites 00288 CALL vsite_force_control(force_env,error) 00289 ! External Potential 00290 CALL add_external_potential(force_env, error=error) 00291 ! Rescale forces if requested 00292 CALL rescale_forces(force_env, error=error) 00293 END IF 00294 IF (ASSOCIATED(force_env%cpot_env)) THEN 00295 CALL force_env_get(force_env, potential_energy=e_pot, error=error) 00296 nat=force_env_get_natom(force_env,error=error) 00297 ALLOCATE(pos(3*nat),stat=stat) 00298 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00299 CALL force_env_get_pos(force_env, pos, 3*nat, error=error) 00300 DEALLOCATE(pos,stat=stat) 00301 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00302 CALL cpot_calc(force_env%cpot_env, e_pot,pos,error) 00303 END IF 00304 00305 ! Print always Energy in the same format for all methods 00306 output_unit = cp_print_key_unit_nr(logger,force_env%force_env_section,"PRINT%PROGRAM_RUN_INFO",& 00307 extension=".Log",error=error) 00308 IF (output_unit > 0) THEN 00309 CALL force_env_get(force_env, potential_energy=e_pot, error=error) 00310 WRITE(output_unit,'(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) energy (a.u.): ",T55,F26.15,/)')& 00311 ADJUSTR(TRIM(use_prog_name(force_env%in_use))),e_pot 00312 END IF 00313 CALL cp_print_key_finished_output(output_unit,logger,force_env%force_env_section,& 00314 "PRINT%PROGRAM_RUN_INFO",error=error) 00315 00316 ! Print Forces, if requested 00317 print_forces = cp_print_key_unit_nr(logger,force_env%force_env_section,"PRINT%FORCES",& 00318 extension=".xyz",error=error) 00319 IF ((print_forces > 0).AND.calculate_forces) THEN 00320 CALL force_env_get(force_env,subsys=subsys,error=error) 00321 CALL cp_subsys_get(subsys,& 00322 core_particles=core_particles,& 00323 particles=particles,& 00324 shell_particles=shell_particles,& 00325 error=error) 00326 ! Variable precision output of the forces 00327 CALL section_vals_val_get(force_env%force_env_section,"PRINT%FORCES%NDIGITS",& 00328 i_val=ndigits,error=error) 00329 CALL write_forces(particles,print_forces,"ATOMIC",ndigits,total_force,error=error) 00330 grand_total_force(:) = total_force(:) 00331 IF (ASSOCIATED(core_particles)) THEN 00332 CALL write_forces(core_particles,print_forces,"CORE",ndigits,total_force,error=error) 00333 grand_total_force(:) = grand_total_force(:) + total_force(:) 00334 END IF 00335 IF (ASSOCIATED(shell_particles)) THEN 00336 CALL write_forces(shell_particles,print_forces,"SHELL",ndigits,total_force,& 00337 grand_total_force,error=error) 00338 END IF 00339 END IF 00340 CALL cp_print_key_finished_output(print_forces,logger,force_env%force_env_section,"PRINT%FORCES",& 00341 error=error) 00342 00343 ! Store results 00344 IF(ASSOCIATED(results))THEN 00345 CALL cp_result_copy(results_in=results,results_out=force_env%results,error=error) 00346 END IF 00347 END IF 00348 00349 ! Write stress tensor 00350 IF (virial%pv_availability) THEN 00351 ! If the virial is defined but we are not computing forces let's zero the 00352 ! virial for consistency 00353 IF (calculate_forces) THEN 00354 output_unit = cp_print_key_unit_nr(logger,force_env%force_env_section,"PRINT%STRESS_TENSOR",& 00355 extension=".stress_tensor",error=error) 00356 IF (output_unit > 0) THEN 00357 CALL section_vals_val_get(force_env%force_env_section,"PRINT%STRESS_TENSOR%NDIGITS",& 00358 i_val=ndigits,error=error) 00359 CALL force_env_get(force_env,cell=cell,error=error) 00360 CALL write_stress_tensor(virial%pv_virial,output_unit,cell,ndigits,virial%pv_numer,& 00361 error=error) 00362 END IF 00363 CALL cp_print_key_finished_output(output_unit,logger,force_env%force_env_section,& 00364 "PRINT%STRESS_TENSOR",error=error) 00365 ELSE 00366 CALL zero_virial(virial,reset=.FALSE.) 00367 END IF 00368 END IF 00369 00370 ! Atomic properties 00371 output_unit = cp_print_key_unit_nr(logger,force_env%force_env_section,"PRINT%PROGRAM_RUN_INFO",& 00372 extension=".Log",error=error) 00373 IF (atprop_env%energy) THEN 00374 CALL mp_sum(atprop_env%atener, force_env%para_env%group) 00375 CALL force_env_get(force_env, potential_energy=e_pot, error=error) 00376 IF (output_unit > 0) THEN 00377 IF (logger%iter_info%print_level > low_print_level) THEN 00378 WRITE (UNIT=output_unit,FMT="(/,T6,A,T15,A)") "Atom","Potential energy" 00379 WRITE (UNIT=output_unit,FMT="(T2,I8,1X,F20.10)")& 00380 (i,atprop_env%atener(i),i=1,SIZE(atprop_env%atener)) 00381 END IF 00382 sum_energy = SUM(atprop_env%atener(:)) 00383 checksum = ABS(e_pot - sum_energy) 00384 WRITE (UNIT=output_unit,FMT="(/,(T2,A,1X,F25.13))")& 00385 "Potential energy (Atomic):",sum_energy,& 00386 "Potential energy (Total) :",e_pot,& 00387 "Difference :",checksum 00388 CPPostcondition((checksum < ateps*ABS(e_pot)),cp_fatal_level,routineP,error,failure) 00389 END IF 00390 END IF 00391 IF (atprop_env%stress) THEN 00392 CALL mp_sum(atprop_env%atstress,force_env%para_env%group) 00393 IF (output_unit > 0) THEN 00394 IF (logger%iter_info%print_level > low_print_level) THEN 00395 DO i=1,SIZE(atprop_env%atstress,3) 00396 WRITE (UNIT=output_unit,FMT="(/,T2,I0,T11,A1,2(14X,A1))") i,"X","Y","Z" 00397 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "X",(atprop_env%atstress(1,j,i),j=1,3) 00398 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Y",(atprop_env%atstress(2,j,i),j=1,3) 00399 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Z",(atprop_env%atstress(3,j,i),j=1,3) 00400 WRITE (UNIT=output_unit,FMT="(T2,A,F15.8)") "1/3 Trace(Atomic stress tensor):",& 00401 (atprop_env%atstress(1,1,i) + atprop_env%atstress(2,2,i) + atprop_env%atstress(3,3,i))/3.0_dp 00402 END DO 00403 END IF 00404 atomic_stress_tensor(:,:) = 0.0_dp 00405 DO i=1,3 00406 atomic_stress_tensor(i,i) = SUM(atprop_env%atstress(i,i,:)) 00407 DO j=i+1,3 00408 atomic_stress_tensor(i,j) = SUM(atprop_env%atstress(i,j,:)) 00409 atomic_stress_tensor(j,i) = atomic_stress_tensor(i,j) 00410 END DO 00411 END DO 00412 WRITE (UNIT=output_unit,FMT="(/,T2,A,T11,A1,2(14X,A1))") "Atomic","X","Y","Z" 00413 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "X",(atomic_stress_tensor(1,i),i=1,3) 00414 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Y",(atomic_stress_tensor(2,i),i=1,3) 00415 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Z",(atomic_stress_tensor(3,i),i=1,3) 00416 WRITE (UNIT=output_unit,FMT="(T2,A,F15.8)") "1/3 Trace(Atomic stress tensor):",& 00417 (atomic_stress_tensor(1,1) + atomic_stress_tensor(2,2) + atomic_stress_tensor(3,3))/3.0_dp 00418 sum_stress_tensor = SUM(atomic_stress_tensor(:,:)) 00419 IF (virial%pv_availability.AND.calculate_forces) THEN 00420 WRITE (UNIT=output_unit,FMT="(/,T2,A,T11,A1,2(14X,A1))") "Total","X","Y","Z" 00421 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "X",(virial%pv_virial(1,i),i=1,3) 00422 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Y",(virial%pv_virial(2,i),i=1,3) 00423 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Z",(virial%pv_virial(3,i),i=1,3) 00424 WRITE (UNIT=output_unit,FMT="(T2,A,F15.8)") "1/3 Trace(Total stress tensor): ",& 00425 (virial%pv_virial(1,1) + virial%pv_virial(2,2) + virial%pv_virial(3,3))/3.0_dp 00426 sum_pv_virial = SUM(virial%pv_virial(:,:)) 00427 diff_stress_tensor(:,:) = ABS(virial%pv_virial(:,:) - atomic_stress_tensor(:,:)) 00428 WRITE (UNIT=output_unit,FMT="(/,T2,A,T11,A1,2(14X,A1))") "Diff","X","Y","Z" 00429 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "X",(diff_stress_tensor(1,i),i=1,3) 00430 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Y",(diff_stress_tensor(2,i),i=1,3) 00431 WRITE (UNIT=output_unit,FMT="(A3,3F15.8)") "Z",(diff_stress_tensor(3,i),i=1,3) 00432 WRITE (UNIT=output_unit,FMT="(T2,A,F15.8)") "1/3 Trace(Diff) : ",& 00433 (diff_stress_tensor(1,1) + diff_stress_tensor(2,2) + diff_stress_tensor(3,3))/3.0_dp 00434 checksum = SUM(diff_stress_tensor(:,:)) 00435 WRITE (UNIT=output_unit,FMT="(/,(T2,A,1X,F25.13))")& 00436 "Checksum stress (Atomic) :",sum_stress_tensor,& 00437 "Checksum stress (Total) :",sum_pv_virial,& 00438 "Difference :",checksum 00439 CPPostcondition((checksum < ateps),cp_fatal_level,routineP,error,failure) 00440 END IF 00441 END IF 00442 END IF 00443 00444 END SUBROUTINE force_env_calc_energy_force 00445 00446 ! ***************************************************************************** 00454 SUBROUTINE force_env_calc_num_pressure(force_env, dx, error) 00455 00456 TYPE(force_env_type), POINTER :: force_env 00457 REAL(KIND=dp), INTENT(IN), OPTIONAL :: dx 00458 TYPE(cp_error_type), INTENT(inout) :: error 00459 00460 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_calc_num_pressure', 00461 routineP = moduleN//':'//routineN 00462 REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp 00463 00464 INTEGER :: i, ip, iq, istat, j, k, 00465 natom, ncore, nshell, 00466 output_unit 00467 LOGICAL :: failure 00468 REAL(KIND=dp) :: dx_w 00469 REAL(KIND=dp), DIMENSION(2) :: numer_energy 00470 REAL(KIND=dp), DIMENSION(3) :: s 00471 REAL(KIND=dp), DIMENSION(3, 3) :: numer_stress = 0.0_dp 00472 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ref_pos_atom, ref_pos_core, 00473 ref_pos_shell 00474 TYPE(cell_type), POINTER :: cell, cell_local 00475 TYPE(cp_logger_type), POINTER :: logger 00476 TYPE(cp_subsys_type), POINTER :: subsys 00477 TYPE(global_environment_type), POINTER :: globenv 00478 TYPE(particle_list_type), POINTER :: core_particles, particles, 00479 shell_particles 00480 TYPE(virial_type), POINTER :: virial 00481 00482 NULLIFY (cell_local) 00483 NULLIFY (core_particles) 00484 NULLIFY (particles) 00485 NULLIFY (shell_particles) 00486 NULLIFY (ref_pos_atom) 00487 NULLIFY (ref_pos_core) 00488 NULLIFY (ref_pos_shell) 00489 natom = 0 00490 ncore = 0 00491 nshell = 0 00492 00493 failure = .FALSE. 00494 logger => cp_error_get_logger(error) 00495 00496 IF (.NOT.failure) THEN 00497 00498 dx_w = default_dx 00499 IF (PRESENT(dx)) dx_w = dx 00500 CALL force_env_get(force_env,subsys=subsys,globenv=globenv,error=error) 00501 CALL cp_subsys_get(subsys,& 00502 core_particles=core_particles,& 00503 particles=particles,& 00504 shell_particles=shell_particles,& 00505 error=error) 00506 output_unit = cp_print_key_unit_nr(logger,force_env%force_env_section,"PRINT%STRESS_TENSOR",& 00507 extension=".stress_tensor",error=error) 00508 IF (output_unit > 0) THEN 00509 WRITE (output_unit,'(/A,A/)') ' **************************** ', & 00510 'NUMERICAL STRESS ********************************' 00511 END IF 00512 00513 ! Save all original particle positions 00514 natom = particles%n_els 00515 ALLOCATE (ref_pos_atom(natom,3),STAT=istat) 00516 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00517 DO i=1,natom 00518 ref_pos_atom(i,:) = particles%els(i)%r 00519 END DO 00520 IF (ASSOCIATED(core_particles)) THEN 00521 ncore = core_particles%n_els 00522 ALLOCATE (ref_pos_core(ncore,3),STAT=istat) 00523 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00524 DO i=1,ncore 00525 ref_pos_core(i,:) = core_particles%els(i)%r 00526 END DO 00527 END IF 00528 IF (ASSOCIATED(shell_particles)) THEN 00529 nshell = shell_particles%n_els 00530 ALLOCATE (ref_pos_shell(nshell,3),STAT=istat) 00531 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00532 DO i=1,nshell 00533 ref_pos_shell(i,:) = shell_particles%els(i)%r 00534 END DO 00535 END IF 00536 CALL force_env_get(force_env,cell=cell,virial=virial,error=error) 00537 CALL cell_create(cell_local,error=error) 00538 CALL cell_clone(cell,cell_local,error=error) 00539 ! First change box 00540 DO ip=1,3 00541 DO iq=1,3 00542 IF (virial%pv_diagonal.AND.(ip /= iq)) CYCLE 00543 DO k=1,2 00544 cell%hmat(ip,iq) = cell_local%hmat(ip,iq) - (-1.0_dp)**k*dx_w 00545 CALL init_cell(cell) 00546 ! Scale positions 00547 DO i=1,natom 00548 CALL real_to_scaled(s,ref_pos_atom(i,1:3),cell_local) 00549 CALL scaled_to_real(particles%els(i)%r,s,cell) 00550 END DO 00551 DO i=1,ncore 00552 CALL real_to_scaled(s,ref_pos_core(i,1:3),cell_local) 00553 CALL scaled_to_real(core_particles%els(i)%r,s,cell) 00554 END DO 00555 DO i=1,nshell 00556 CALL real_to_scaled(s,ref_pos_shell(i,1:3),cell_local) 00557 CALL scaled_to_real(shell_particles%els(i)%r,s,cell) 00558 END DO 00559 ! Since the box has changed, rebuild grids, i.e. pw_env and ks_env 00560 CALL force_env_set_cell(force_env,cell=cell,error=error) 00561 ! Compute energies 00562 CALL force_env_calc_energy_force(force_env,calc_force=.FALSE.,& 00563 consistent_energies=.TRUE.,error=error) 00564 CALL force_env_get(force_env,potential_energy=numer_energy(k),error=error) 00565 ! Reset cell 00566 cell%hmat(ip,iq) = cell_local%hmat(ip,iq) 00567 END DO 00568 CALL init_cell(cell) 00569 numer_stress(ip,iq) = (numer_energy(1) - numer_energy(2))/(2.0_dp*dx_w) 00570 END DO 00571 END DO 00572 00573 ! Reset positions and rebuild original environment 00574 DO i=1,natom 00575 particles%els(i)%r = ref_pos_atom(i,:) 00576 END DO 00577 DO i=1,ncore 00578 core_particles%els(i)%r = ref_pos_core(i,:) 00579 END DO 00580 DO i=1,nshell 00581 shell_particles%els(i)%r = ref_pos_shell(i,:) 00582 END DO 00583 CALL force_env_set_cell(force_env,cell=cell,error=error) 00584 CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.,& 00585 consistent_energies=.TRUE.,& 00586 error=error) 00587 00588 ! Computing pv_test 00589 virial%pv_virial = 0.0_dp 00590 DO i=1,3 00591 DO j=1,3 00592 DO k=1,3 00593 virial%pv_virial(i,j) = virial%pv_virial(i,j) -& 00594 0.5_dp*(numer_stress(i,k)*cell_local%hmat(j,k) +& 00595 numer_stress(j,k)*cell_local%hmat(i,k)) 00596 END DO 00597 END DO 00598 END DO 00599 IF (globenv%run_type_id == debug_run) THEN 00600 IF (output_unit > 0) THEN 00601 WRITE (UNIT=output_unit,FMT="((T2,A))") "Numerical pv_virial" 00602 WRITE (UNIT=output_unit,FMT="((T3,3F16.10))") (virial%pv_virial(i,1:3),i=1,3) 00603 END IF 00604 END IF 00605 00606 IF (output_unit > 0) THEN 00607 WRITE (output_unit,'(/,A,/)') ' **************************** '//& 00608 'NUMERICAL STRESS END *****************************' 00609 END IF 00610 00611 CALL cp_print_key_finished_output(output_unit,logger,force_env%force_env_section,& 00612 "PRINT%STRESS_TENSOR",error=error) 00613 00614 ! Release storage 00615 IF (ASSOCIATED(ref_pos_atom)) THEN 00616 DEALLOCATE (ref_pos_atom,STAT=istat) 00617 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00618 END IF 00619 IF (ASSOCIATED(ref_pos_core)) THEN 00620 DEALLOCATE (ref_pos_core,STAT=istat) 00621 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00622 END IF 00623 IF (ASSOCIATED(ref_pos_shell)) THEN 00624 DEALLOCATE (ref_pos_shell,STAT=istat) 00625 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00626 END IF 00627 IF (ASSOCIATED(cell_local)) CALL cell_release(cell_local,error=error) 00628 00629 END IF 00630 00631 END SUBROUTINE force_env_calc_num_pressure 00632 00633 ! ***************************************************************************** 00644 SUBROUTINE force_env_create(force_env,root_section,para_env,globenv,fist_env,& 00645 qs_env,meta_env,sub_force_env,qmmm_env,eip_env,ep_env,force_env_section,& 00646 mixed_env,error) 00647 00648 TYPE(force_env_type), POINTER :: force_env 00649 TYPE(section_vals_type), POINTER :: root_section 00650 TYPE(cp_para_env_type), POINTER :: para_env 00651 TYPE(global_environment_type), POINTER :: globenv 00652 TYPE(fist_environment_type), OPTIONAL, 00653 POINTER :: fist_env 00654 TYPE(qs_environment_type), OPTIONAL, 00655 POINTER :: qs_env 00656 TYPE(meta_env_type), OPTIONAL, POINTER :: meta_env 00657 TYPE(force_env_p_type), DIMENSION(:), 00658 OPTIONAL, POINTER :: sub_force_env 00659 TYPE(qmmm_env_qm_type), OPTIONAL, 00660 POINTER :: qmmm_env 00661 TYPE(eip_environment_type), OPTIONAL, 00662 POINTER :: eip_env 00663 TYPE(ep_env_type), OPTIONAL, POINTER :: ep_env 00664 TYPE(section_vals_type), POINTER :: force_env_section 00665 TYPE(mixed_environment_type), OPTIONAL, 00666 POINTER :: mixed_env 00667 TYPE(cp_error_type), INTENT(inout) :: error 00668 00669 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_create', 00670 routineP = moduleN//':'//routineN 00671 00672 INTEGER :: stat, stress_tensor 00673 LOGICAL :: atomic_energy, atomic_stress, 00674 failure, pv_availability, 00675 pv_diagonal, pv_numerical 00676 TYPE(cp_subsys_type), POINTER :: subsys 00677 00678 failure=.FALSE. 00679 00680 ALLOCATE ( force_env, stat=stat ) 00681 CPPostconditionNoFail(stat==0,cp_fatal_level,routineP,error) 00682 IF (.NOT. failure) THEN 00683 NULLIFY ( force_env%subsys, force_env%fist_env, & 00684 force_env%qs_env, & 00685 force_env%para_env, force_env%globenv, & 00686 force_env%meta_env, force_env%sub_force_env, & 00687 force_env%qmmm_env, force_env%ep_env, force_env%fp_env, & 00688 force_env%force_env_section, force_env%eip_env,force_env%mixed_env,& 00689 force_env%root_section, force_env%cpot_env, force_env%atprop_env, force_env%results) 00690 last_force_env_id=last_force_env_id+1 00691 force_env%id_nr=last_force_env_id 00692 force_env%ref_count=1 00693 force_env%in_use=0 00694 force_env%additional_potential=0.0_dp 00695 00696 force_env%globenv => globenv 00697 CALL globenv_retain(force_env%globenv,error=error) 00698 00699 force_env%root_section => root_section 00700 CALL section_vals_retain(root_section,error=error) 00701 00702 force_env%para_env=>para_env 00703 CALL cp_para_env_retain(force_env%para_env, error=error) 00704 00705 CALL section_vals_retain(force_env_section,error=error) 00706 force_env%force_env_section => force_env_section 00707 00708 ! Should we compute the virial? 00709 CALL section_vals_val_get(force_env_section,"STRESS_TENSOR",i_val=stress_tensor,error=error) 00710 SELECT CASE(stress_tensor) 00711 CASE(do_stress_none) 00712 pv_availability=.FALSE. 00713 pv_numerical=.FALSE. 00714 pv_diagonal=.FALSE. 00715 CASE(do_stress_analytical) 00716 pv_availability=.TRUE. 00717 pv_numerical=.FALSE. 00718 pv_diagonal=.FALSE. 00719 CASE(do_stress_numerical) 00720 pv_availability=.TRUE. 00721 pv_numerical=.TRUE. 00722 pv_diagonal=.FALSE. 00723 CASE(do_stress_diagonal_anal) 00724 pv_availability=.TRUE. 00725 pv_numerical=.FALSE. 00726 pv_diagonal=.TRUE. 00727 CASE(do_stress_diagonal_numer) 00728 pv_availability=.TRUE. 00729 pv_numerical=.TRUE. 00730 pv_diagonal=.TRUE. 00731 END SELECT 00732 00733 ! Should we compute atomic properties? 00734 CALL atprop_create(force_env%atprop_env,error) 00735 CALL section_vals_val_get(force_env_section,"PROPERTIES%ATOMIC%ENERGY",l_val=atomic_energy,error=error) 00736 force_env%atprop_env%energy = atomic_energy 00737 CALL section_vals_val_get(force_env_section,"PROPERTIES%ATOMIC%PRESSURE",l_val=atomic_stress,error=error) 00738 IF (atomic_stress) THEN 00739 CPPrecondition(pv_availability,cp_failure_level,routineP,error,failure) 00740 CPPrecondition(.NOT.pv_numerical,cp_failure_level,routineP,error,failure) 00741 END IF 00742 force_env%atprop_env%stress = atomic_stress 00743 00744 IF (PRESENT(fist_env)) THEN 00745 IF (ASSOCIATED(fist_env)) THEN 00746 CPPrecondition(force_env%in_use==0,cp_failure_level,routineP,error,failure) 00747 force_env%in_use=use_fist_force 00748 force_env%fist_env => fist_env 00749 CALL fist_env_retain(fist_env,error=error) 00750 ! Virial controlled through the external request 00751 CALL virial_create(force_env%virial,error=error) 00752 CALL virial_set(force_env%virial,& 00753 pv_availability=pv_availability,& 00754 pv_numer=pv_numerical,& 00755 pv_diagonal=pv_diagonal) 00756 END IF 00757 END IF 00758 IF (PRESENT(eip_env)) THEN 00759 IF (ASSOCIATED(eip_env)) THEN 00760 CPPrecondition(force_env%in_use==0, cp_failure_level, routineP, error, failure) 00761 force_env%in_use = use_eip_force 00762 force_env%eip_env => eip_env 00763 CALL eip_env_retain(eip_env, error=error) 00764 ! Virial not present for EIP 00765 CALL virial_create(force_env%virial, error=error) 00766 eip_env%virial => force_env%virial 00767 CALL virial_retain(eip_env%virial,error=error) 00768 END IF 00769 END IF 00770 IF (PRESENT(qs_env)) THEN 00771 IF (ASSOCIATED(qs_env)) THEN 00772 CPPrecondition(force_env%in_use==0,cp_failure_level,routineP,error,failure) 00773 force_env%in_use=use_qs_force 00774 force_env%qs_env => qs_env 00775 CALL qs_env_retain(qs_env,error=error) 00776 CALL virial_create(force_env%virial, error=error) 00777 ! Virial controlled through the external request 00778 CALL virial_set(virial=force_env%virial,& 00779 pv_availability=pv_availability,& 00780 pv_numer=pv_numerical,& 00781 pv_diagonal=pv_diagonal) 00782 qs_env%virial => force_env%virial 00783 CALL virial_retain(qs_env%virial,error=error) 00784 END IF 00785 END IF 00786 IF (PRESENT(qmmm_env)) THEN 00787 CPPrecondition(PRESENT(sub_force_env),cp_failure_level,routineP,error,failure) 00788 force_env%in_use=use_qmmm 00789 force_env%qmmm_env => qmmm_env 00790 CALL qmmm_env_qm_retain(qmmm_env,error=error) 00791 force_env%virial => sub_force_env(primary_subsys)%force_env%virial 00792 CALL virial_retain(force_env%virial,error=error) 00793 ! Virial controlled through the external request 00794 CALL virial_set ( virial=force_env%virial,& 00795 pv_availability=pv_availability,& 00796 pv_numer=pv_numerical,& 00797 pv_diagonal=pv_diagonal) 00798 END IF 00799 IF (PRESENT(mixed_env)) THEN 00800 CPPrecondition(force_env%in_use==0, cp_failure_level, routineP, error, failure) 00801 force_env%in_use=use_mixed_force 00802 force_env%mixed_env => mixed_env 00803 CALL mixed_env_retain ( mixed_env, error = error ) 00804 ! This is necessary as long as there are methods not implementing the virial 00805 CALL virial_create ( force_env % virial, error=error) 00806 CALL virial_set ( virial=force_env%virial,& 00807 pv_availability=pv_availability,& 00808 pv_numer=pv_numerical,& 00809 pv_diagonal=pv_diagonal) 00810 END IF 00811 IF (PRESENT(ep_env)) THEN 00812 IF (ASSOCIATED(ep_env)) THEN 00813 CPPrecondition(force_env%in_use==0,cp_failure_level,routineP,error,failure) 00814 force_env%in_use=use_ep_force 00815 force_env%ep_env => ep_env 00816 CALL ep_env_retain(ep_env,error=error) 00817 ! Virial not present for EP 00818 CALL virial_create ( force_env%virial, error=error) 00819 END IF 00820 END IF 00821 CPPostcondition(force_env%in_use/=0,cp_failure_level,routineP,error,failure) 00822 00823 IF (PRESENT(sub_force_env)) THEN 00824 force_env%sub_force_env => sub_force_env 00825 END IF 00826 00827 IF (PRESENT(meta_env)) THEN 00828 force_env%meta_env => meta_env 00829 CALL meta_env_retain(meta_env,error=error) 00830 ELSE 00831 NULLIFY(force_env%meta_env) 00832 END IF 00833 CALL cp_result_create(results=force_env%results,error=error) 00834 SELECT CASE(force_env%in_use) 00835 CASE(use_fist_force) 00836 CALL fist_env_get (force_env%fist_env, subsys=force_env%subsys, error=error) 00837 CALL cp_subsys_retain (force_env%subsys, error=error) 00838 CASE(use_qs_force) 00839 CALL get_qs_env(force_env%qs_env, subsys=force_env%subsys,error=error) 00840 CALL cp_subsys_retain(force_env%subsys,error=error) 00841 CASE(use_ep_force) 00842 CALL get_qs_env(force_env%ep_env%main_qs_env,subsys=force_env%subsys,error=error) 00843 CALL cp_subsys_retain(force_env%subsys,error=error) 00844 CASE(use_eip_force) 00845 CALL eip_env_get(force_env%eip_env,subsys=force_env%subsys,error=error) 00846 CALL cp_subsys_retain(force_env%subsys, error=error) 00847 CASE(use_qmmm) 00848 subsys => force_env%sub_force_env(primary_subsys)%force_env%subsys 00849 force_env%subsys => subsys 00850 CALL cp_subsys_retain(subsys,error=error) 00851 CASE(use_mixed_force) 00852 CALL get_mixed_env (force_env%mixed_env,subsys=force_env%subsys,error=error) 00853 CALL cp_subsys_retain (force_env%subsys, error=error) 00854 CASE default 00855 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00856 END SELECT 00857 00858 END IF 00859 00860 END SUBROUTINE force_env_create 00861 00862 ! ***************************************************************************** 00870 SUBROUTINE ep_create_force_env(force_env, root_section, para_env, globenv,& 00871 force_env_section, error) 00872 00873 TYPE(force_env_type), POINTER :: force_env 00874 TYPE(section_vals_type), POINTER :: root_section 00875 TYPE(cp_para_env_type), POINTER :: para_env 00876 TYPE(global_environment_type), POINTER :: globenv 00877 TYPE(section_vals_type), POINTER :: force_env_section 00878 TYPE(cp_error_type), INTENT(inout) :: error 00879 00880 CHARACTER(len=*), PARAMETER :: routineN = 'ep_create_force_env', 00881 routineP = moduleN//':'//routineN 00882 00883 LOGICAL :: failure 00884 TYPE(ep_env_type), POINTER :: ep_env 00885 00886 failure=.FALSE. 00887 00888 IF (.NOT. failure) THEN 00889 NULLIFY(ep_env) 00890 CALL ep_env_create(ep_env, root_section, para_env, globenv=globenv,& 00891 error=error) 00892 CALL force_env_create(force_env,root_section,para_env,globenv=globenv,ep_env=ep_env,& 00893 force_env_section = force_env_section, error=error) 00894 CALL ep_env_release(ep_env,error=error) 00895 END IF 00896 END SUBROUTINE ep_create_force_env 00897 00898 ! ***************************************************************************** 00907 RECURSIVE SUBROUTINE qmmm_energy_and_forces(force_env,calc_force,require_consistent_energy_force,error) 00908 TYPE(force_env_type), POINTER :: force_env 00909 LOGICAL, INTENT(IN), OPTIONAL :: calc_force, 00910 require_consistent_energy_force 00911 TYPE(cp_error_type), INTENT(inout) :: error 00912 00913 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces', 00914 routineP = moduleN//':'//routineN 00915 00916 CHARACTER(DEFAULT_STRING_LENGTH) :: restart_filename, 00917 restart_filename_suffix, 00918 restart_history_filename 00919 INTEGER :: ip, isubf, 00920 mom_conserv_min_label, 00921 mom_conserv_n, 00922 mom_conserv_region, 00923 mom_conserv_type 00924 INTEGER, POINTER :: cur_indices(:), cur_labels(:) 00925 LOGICAL :: do_require_consistent_energy_force, failure 00926 REAL(dp) :: delta_a(3), delta_f(3), 00927 mom_conserv_mass, total_f(3) 00928 TYPE(cp_subsys_type), POINTER :: subsys_primary, 00929 subsys_qmmm_core, 00930 subsys_qmmm_extended 00931 TYPE(particle_type), DIMENSION(:), 00932 POINTER :: particles_primary, 00933 particles_qmmm_core, 00934 particles_qmmm_extended 00935 TYPE(section_vals_type), POINTER :: force_env_section 00936 00937 do_require_consistent_energy_force = .TRUE. 00938 IF (PRESENT(require_consistent_energy_force)) do_require_consistent_energy_force = require_consistent_energy_force 00939 00940 IF (do_require_consistent_energy_force) THEN 00941 CALL cp_assert(.NOT. qmmm_force_mixing_active(force_env, error),& 00942 cp_failure_level,cp_assertion_failed,& 00943 routineP,"qmmm_energy_and_forces got require_consistent_energy_force but force mixing is active. "//& 00944 CPSourceFileRef,& 00945 error,failure) 00946 ENDIF 00947 00948 ! Possibly translate the system 00949 00950 CALL apply_qmmm_translate(force_env, error) 00951 00952 DO isubf=1, SIZE(force_env%sub_force_env) 00953 IF (SIZE(force_env%sub_force_env) > 1) THEN ! rewrite restart file name 00954 CALL force_env_get(force_env, force_env_section=force_env_section,error=error) 00955 CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", & 00956 c_val=restart_filename, error=error) 00957 CALL section_vals_val_get(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", & 00958 c_val=restart_history_filename, error=error) 00959 00960 WRITE (unit=restart_filename_suffix, fmt=*) isubf 00961 restart_filename_suffix= ADJUSTL(restart_filename_suffix) 00962 00963 CALL force_env_get(force_env%sub_force_env(isubf)%force_env, force_env_section=force_env_section,error=error) 00964 CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", & 00965 c_val=TRIM(restart_filename)//"-SubForceEnv-"//TRIM(restart_filename_suffix), error=error) 00966 CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", & 00967 c_val=TRIM(restart_history_filename)//"-SubForceEnv-"//TRIM(restart_filename_suffix), error=error) 00968 ENDIF 00969 00970 CALL qmmm_energy_and_forces_low(force_env%sub_force_env(isubf)%force_env,& 00971 calc_force,error) 00972 00973 IF (SIZE(force_env%sub_force_env) > 1) THEN ! undo rewrite of restart file name 00974 CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART%FILENAME", & 00975 c_val=TRIM(restart_filename), error=error) 00976 CALL section_vals_val_set(force_env_section, "DFT%SCF%PRINT%RESTART_HISTORY%FILENAME", & 00977 c_val=TRIM(restart_history_filename), error=error) 00978 ENDIF 00979 END DO 00980 00981 IF (SIZE(force_env%sub_force_env) == 1) THEN 00982 ! see force_env_get for potential energy to copy 00983 ! copy forces, energy, virial up 00984 CONTINUE 00985 ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN 00986 ! force mixing happens here 00987 ! energy will just inherit from extended sub_force_env in force_env_get() 00988 00989 ! get forces from subsys of each sub force env 00990 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,& 00991 subsys=subsys_qmmm_core,force_env_section=force_env_section,error=error) 00992 CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,& 00993 subsys=subsys_qmmm_extended,error=error) 00994 00995 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%INDICES",i_vals=cur_indices,error=error) 00996 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%RESTART_INFO%LABELS",i_vals=cur_labels,error=error) 00997 00998 particles_qmmm_extended => subsys_qmmm_extended%particles%els 00999 particles_qmmm_core => subsys_qmmm_core%particles%els 01000 DO ip=1,SIZE(cur_indices) 01001 IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom 01002 ! copy (QM) force from extended calculation 01003 particles_qmmm_core(cur_indices(ip))%f=particles_qmmm_extended(cur_indices(ip))%f 01004 ENDIF 01005 END DO 01006 01007 ! zero momentum 01008 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_TYPE",& 01009 i_val=mom_conserv_type,error=error) 01010 IF (mom_conserv_type /= do_fm_mom_conserv_none) THEN 01011 CALL section_vals_val_get(force_env_section,"QMMM%FORCE_MIXING%MOMENTUM_CONSERVATION_REGION",& 01012 i_val=mom_conserv_region,error=error) 01013 01014 IF (mom_conserv_region == do_fm_mom_conserv_core) THEN 01015 mom_conserv_min_label = force_mixing_label_QM_core 01016 ELSEIF (mom_conserv_region == do_fm_mom_conserv_QM) THEN 01017 mom_conserv_min_label = force_mixing_label_QM_dynamics 01018 ELSEIF (mom_conserv_region == do_fm_mom_conserv_buffer) THEN 01019 mom_conserv_min_label = force_mixing_label_buffer 01020 ELSE 01021 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01022 routineP,"Got unknown MOMENTUM_CONSERVATION_REGION (not CORE, QM, or BUFFER) !"//& 01023 CPSourceFileRef,& 01024 error,failure) 01025 ENDIF 01026 01027 total_f = 0.0_dp 01028 DO ip=1, SIZE(particles_qmmm_core) 01029 total_f(1:3) = total_f(1:3) + particles_qmmm_core(ip)%f(1:3) 01030 END DO 01031 IF (mom_conserv_type == do_fm_mom_conserv_equal_f) THEN 01032 mom_conserv_n = COUNT(cur_labels >= mom_conserv_region) 01033 delta_f = total_f/mom_conserv_n 01034 DO ip=1, SIZE(cur_indices) 01035 IF (cur_labels(ip) >= mom_conserv_region) THEN 01036 particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - delta_f 01037 ENDIF 01038 END DO 01039 ELSE IF (mom_conserv_type == do_fm_mom_conserv_equal_a) THEN 01040 mom_conserv_mass = 0.0_dp 01041 DO ip=1, SIZE(cur_indices) 01042 IF (cur_labels(ip) >= mom_conserv_region) & 01043 mom_conserv_mass = mom_conserv_mass + particles_qmmm_core(cur_indices(ip))%atomic_kind%mass 01044 END DO 01045 delta_a = total_f/mom_conserv_mass 01046 DO ip=1, SIZE(cur_indices) 01047 IF (cur_labels(ip) >= mom_conserv_region) THEN 01048 particles_qmmm_core(cur_indices(ip))%f = particles_qmmm_core(cur_indices(ip))%f - & 01049 particles_qmmm_core(cur_indices(ip))%atomic_kind%mass * delta_a 01050 ENDIF 01051 END DO 01052 ENDIF 01053 ENDIF 01054 01055 IF (force_mixing_core_subsys /= primary_subsys) THEN 01056 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env,& 01057 subsys=subsys_primary,force_env_section=force_env_section,error=error) 01058 particles_primary => subsys_primary%particles%els 01059 DO ip=1,SIZE(particles_qmmm_core) 01060 particles_primary(ip)%f=particles_qmmm_core(ip)%f 01061 END DO 01062 ENDIF 01063 01064 ELSE 01065 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01066 routineP,"size(sub_force_env) neither 1 (conventional QM/MM) nor "//& 01067 "2 (force-mixing)! "//& 01068 CPSourceFileRef,& 01069 error,failure) 01070 ENDIF 01071 01072 END SUBROUTINE qmmm_energy_and_forces 01073 01074 ! ***************************************************************************** 01083 RECURSIVE SUBROUTINE qmmm_energy_and_forces_low(force_env,calc_force,error) 01084 01085 TYPE(force_env_type), POINTER :: force_env 01086 LOGICAL, INTENT(IN), OPTIONAL :: calc_force 01087 TYPE(cp_error_type), INTENT(inout) :: error 01088 01089 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_energy_and_forces_low', 01090 routineP = moduleN//':'//routineN 01091 01092 INTEGER :: ip, output_unit 01093 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 01094 LOGICAL :: calculate_forces, failure, 01095 qmmm_added_chrg, qmmm_link, 01096 qmmm_link_imomm 01097 REAL(KIND=dp) :: energy_mm, energy_qm 01098 REAL(KIND=dp), DIMENSION(3) :: max_coord, min_coord 01099 TYPE(cell_type), POINTER :: mm_cell, qm_cell 01100 TYPE(cp_logger_type), POINTER :: logger 01101 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 01102 TYPE(particle_type), DIMENSION(:), 01103 POINTER :: particles_mm, particles_qm 01104 TYPE(qmmm_links_type), POINTER :: qmmm_links 01105 TYPE(section_vals_type), POINTER :: force_env_section 01106 01107 min_coord = HUGE(0.0_dp) 01108 max_coord = -HUGE(0.0_dp) 01109 failure = .FALSE. 01110 calculate_forces = .TRUE. 01111 qmmm_link = .FALSE. 01112 qmmm_link_imomm = .FALSE. 01113 qmmm_added_chrg = .FALSE. 01114 logger => cp_error_get_logger(error) 01115 IF (PRESENT(calc_force)) calculate_forces = calc_force 01116 NULLIFY(subsys_mm, subsys_qm, qm_atom_index,particles_mm,particles_qm, qm_cell, mm_cell) 01117 NULLIFY(force_env_section) 01118 force_env_section => force_env%sub_force_env(qs_subsys)%force_env%force_env_section 01119 01120 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 01121 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 01122 CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) 01123 CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) 01124 01125 CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,& 01126 cell=mm_cell,subsys=subsys_mm,error=error) 01127 CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& 01128 cell=qm_cell,subsys=subsys_qm,error=error) 01129 qm_atom_index => force_env%qmmm_env%qm_atom_index 01130 qmmm_link = force_env%qmmm_env%qmmm_link 01131 qmmm_links => force_env%qmmm_env%qmmm_links 01132 qmmm_added_chrg = (force_env%qmmm_env%move_mm_charges .OR. force_env%qmmm_env%add_mm_charges) 01133 IF (qmmm_link) THEN 01134 CPPrecondition(ASSOCIATED(qmmm_links),cp_failure_level,routineP,error,failure) 01135 IF (ASSOCIATED(qmmm_links%imomm)) qmmm_link_imomm = (SIZE(qmmm_links%imomm) /= 0) 01136 END IF 01137 CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) 01138 01139 particles_mm => subsys_mm%particles%els 01140 particles_qm => subsys_qm%particles%els 01141 01142 ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom 01143 IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, particles_qm, qm_atom_index, error) 01144 01145 ! If add charges get their position NOW! 01146 IF (qmmm_added_chrg) CALL qmmm_added_chrg_coord(force_env%qmmm_env, particles_mm, error) 01147 01148 ! Initialize ks_qmmm_env 01149 CALL ks_qmmm_env_rebuild(qs_env=force_env%sub_force_env(qs_subsys)%force_env%qs_env,& 01150 qmmm_env=force_env%qmmm_env,error=error) 01151 01152 ! Compute the short range QM/MM Electrostatic Potential 01153 CALL qmmm_el_coupling( qs_env=force_env%sub_force_env(qs_subsys)%force_env%qs_env,& 01154 qmmm_env=force_env%qmmm_env,& 01155 mm_particles=particles_mm,& 01156 mm_cell=mm_cell,& 01157 error=error) 01158 01159 ! Fist 01160 CALL force_env_calc_energy_force(force_env%sub_force_env(fist_subsys)%force_env,& 01161 calc_force=calculate_forces,skip_external_control=.TRUE.,error=error) 01162 01163 ! Print Out information on fist energy calculation... 01164 CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,& 01165 potential_energy=energy_mm,& 01166 error=error) 01167 ! QS 01168 CALL force_env_calc_energy_force(force_env%sub_force_env(qs_subsys)%force_env,& 01169 calc_force=calculate_forces,skip_external_control=.TRUE.,error=error) 01170 ! Print Out information on QS energy calculation... 01171 CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& 01172 potential_energy=energy_qm,& 01173 error=error) 01174 ! QM/MM Interaction Potential forces 01175 CALL qmmm_forces(force_env%sub_force_env(qs_subsys)%force_env%qs_env,& 01176 force_env%qmmm_env,particles_mm,& 01177 mm_cell=mm_cell,& 01178 calc_force=calculate_forces,& 01179 error=error) 01180 ! Forces of quadratic wall on QM atoms 01181 CALL apply_qmmm_walls(force_env,error) 01182 01183 ! Print Out information on QS energy calculation... 01184 CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& 01185 potential_energy=energy_qm,& 01186 error=error) 01187 01188 IF (calculate_forces) THEN 01189 ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom 01190 IF (qmmm_link_imomm) CALL qmmm_link_Imomm_forces(qmmm_links,particles_qm,qm_atom_index,error) 01191 particles_mm => subsys_mm%particles%els 01192 DO ip=1,SIZE(qm_atom_index) 01193 particles_mm(qm_atom_index(ip))%f=particles_mm(qm_atom_index(ip))%f+particles_qm(ip)%f 01194 END DO 01195 ! If add charges get rid of their derivatives right NOW! 01196 IF (qmmm_added_chrg) CALL qmmm_added_chrg_forces(force_env%qmmm_env, particles_mm, error) 01197 END IF 01198 01199 output_unit = cp_print_key_unit_nr(logger,force_env_section,"QMMM%PRINT%DERIVATIVES",& 01200 extension=".Log",error=error) 01201 IF (output_unit>0) THEN 01202 WRITE (unit=output_unit,fmt='(/1X,A,F15.9)')"Energy after QMMM calculation: ",energy_qm 01203 IF (calculate_forces) THEN 01204 WRITE (unit=output_unit,fmt='(/1X,A)')"Derivatives on all atoms after QMMM calculation: " 01205 DO ip=1,SIZE(particles_mm) 01206 WRITE (unit=output_unit,fmt='(1X,3F15.9)') particles_mm(ip)%f 01207 END DO 01208 END IF 01209 END IF 01210 CALL cp_print_key_finished_output(output_unit,logger,force_env_section,& 01211 "QMMM%PRINT%DERIVATIVES",error=error) 01212 END SUBROUTINE qmmm_energy_and_forces_low 01213 01214 ! ***************************************************************************** 01226 SUBROUTINE mixed_energy_forces(force_env, calculate_forces, error) 01227 01228 TYPE(force_env_type), POINTER :: force_env 01229 LOGICAL, INTENT(IN) :: calculate_forces 01230 TYPE(cp_error_type), INTENT(inout) :: error 01231 01232 CHARACTER(len=*), PARAMETER :: routineN = 'mixed_energy_forces', 01233 routineP = moduleN//':'//routineN 01234 01235 CHARACTER(LEN=default_path_length) :: coupling_function 01236 CHARACTER(LEN=default_string_length) :: def_error, description, 01237 this_error 01238 INTEGER :: iforce_eval, iparticle, jparticle, mixing_type, my_group, 01239 natom, nforce_eval, source, stat, unit_nr 01240 INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index 01241 LOGICAL :: failure, virial_consistent 01242 REAL(KIND=dp) :: coupling_parameter, dedf, der_1, der_2, dx, energy, err, 01243 lambda, lerr, restraint_strength, restraint_target, sd 01244 REAL(KIND=dp), DIMENSION(3) :: dip_mix 01245 REAL(KIND=dp), DIMENSION(:), POINTER :: energies 01246 TYPE(cell_type), POINTER :: cell, cell_mix 01247 TYPE(cp_error_type) :: my_error 01248 TYPE(cp_logger_type), POINTER :: logger, my_logger 01249 TYPE(cp_result_p_type), DIMENSION(:), 01250 POINTER :: results 01251 TYPE(cp_result_type), POINTER :: loc_results, results_mix 01252 TYPE(cp_subsys_p_type), DIMENSION(:), 01253 POINTER :: subsystems 01254 TYPE(cp_subsys_type), POINTER :: subsys_mix 01255 TYPE(mixed_energy_type), POINTER :: mixed_energy 01256 TYPE(mixed_force_type), DIMENSION(:), 01257 POINTER :: global_forces 01258 TYPE(particle_list_p_type), 01259 DIMENSION(:), POINTER :: particles 01260 TYPE(particle_list_type), POINTER :: particles_mix 01261 TYPE(section_vals_type), POINTER :: force_env_section, 01262 gen_section, mapping_section, 01263 mixed_section, root_section 01264 TYPE(virial_p_type), DIMENSION(:), 01265 POINTER :: virials 01266 TYPE(virial_type), POINTER :: loc_virial, virial_mix 01267 01268 failure=.FALSE. 01269 logger => cp_error_get_logger(error) 01270 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 01271 ! Get infos about the mixed subsys 01272 CALL force_env_get(force_env=force_env,& 01273 subsys=subsys_mix,& 01274 force_env_section=force_env_section,& 01275 root_section=root_section,& 01276 virial=virial_mix,& 01277 results=results_mix,& 01278 cell=cell_mix,& 01279 error=error) 01280 CALL cp_subsys_get(subsys=subsys_mix,& 01281 particles=particles_mix,& 01282 error=error) 01283 NULLIFY(map_index, glob_natoms, global_forces) 01284 virial_consistent = .TRUE. 01285 nforce_eval = SIZE(force_env%sub_force_env) 01286 mixed_section => section_vals_get_subs_vals(force_env_section,"MIXED",error=error) 01287 mapping_section => section_vals_get_subs_vals(mixed_section,"MAPPING",error=error) 01288 ! Global Info 01289 ALLOCATE(subsystems(nforce_eval), stat=stat) 01290 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01291 ALLOCATE(particles(nforce_eval), stat=stat) 01292 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01293 ! Local Info to sync 01294 ALLOCATE(global_forces(nforce_eval), stat=stat) 01295 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01296 ALLOCATE(energies(nforce_eval), stat=stat) 01297 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01298 ALLOCATE(glob_natoms(nforce_eval), stat=stat) 01299 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01300 ALLOCATE(virials(nforce_eval), stat=stat) 01301 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01302 ALLOCATE(results(nforce_eval), stat=stat) 01303 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01304 energies = 0.0_dp 01305 glob_natoms = 0 01306 DO iforce_eval = 1, nforce_eval 01307 NULLIFY(subsystems(iforce_eval)%subsys, particles(iforce_eval)%list) 01308 NULLIFY(results(iforce_eval)%results, virials(iforce_eval)%virial) 01309 CALL virial_create (virials(iforce_eval)%virial, error) 01310 CALL cp_result_create (results(iforce_eval)%results, error) 01311 IF (.NOT.ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE 01312 ! From this point on the error is the sub_error 01313 my_group = force_env%mixed_env%group_distribution(force_env%para_env%mepos) 01314 my_error = force_env%mixed_env%sub_error(my_group+1) 01315 my_logger => cp_error_get_logger(my_error) 01316 ! Copy iterations info (they are updated only in the main mixed_env) 01317 CALL cp_iteration_info_copy_iter(logger%iter_info, my_logger%iter_info) 01318 01319 ! Get all available subsys 01320 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env,& 01321 subsys=subsystems(iforce_eval)%subsys,cell=cell,error=my_error) 01322 ! Check whether virial can be consistently used.. 01323 IF (virial_mix%pv_availability) THEN 01324 virial_consistent = virial_consistent.AND.compare_cells(cell_mix, cell, my_error) 01325 END IF 01326 ! Get available particles 01327 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys,& 01328 particles=particles(iforce_eval)%list,error=my_error) 01329 01330 ! Get Mapping index array 01331 natom = SIZE(particles(iforce_eval)%list%els) 01332 CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, & 01333 map_index, my_error) 01334 01335 ! Mapping particles from iforce_eval environment to the mixed env 01336 DO iparticle = 1, natom 01337 jparticle = map_index(iparticle) 01338 particles(iforce_eval)%list%els(iparticle)%r= particles_mix%els(jparticle)%r 01339 END DO 01340 01341 ! Calculate energy and forces for each sub_force_env 01342 CALL force_env_calc_energy_force(force_env%sub_force_env(iforce_eval)%force_env,& 01343 calc_force=calculate_forces,& 01344 skip_external_control=.TRUE.,& 01345 error=my_error) 01346 ! Only the rank 0 process collect info for each computation 01347 IF ( force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos==& 01348 force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN 01349 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env,& 01350 potential_energy=energy,& 01351 virial=loc_virial,& 01352 results=loc_results,& 01353 error=my_error) 01354 energies(iforce_eval) = energy 01355 glob_natoms(iforce_eval) = natom 01356 CALL cp_virial(loc_virial, virials(iforce_eval)%virial) 01357 CALL cp_result_copy(loc_results, results(iforce_eval)%results, error) 01358 END IF 01359 ! Deallocate map_index array 01360 IF (ASSOCIATED(map_index)) THEN 01361 DEALLOCATE(map_index, stat=stat) 01362 CPPrecondition(stat==0,cp_failure_level,routineP,my_error,failure) 01363 END IF 01364 CALL cp_error_check(my_error, failure) 01365 END DO 01366 ! Final check on virial 01367 CALL cp_assert(virial_consistent,cp_failure_level,cp_assertion_failed,& 01368 routineP,"Mixed force_eval have different cells definition. Virial cannot be "//& 01369 " defined in a consistent way. Check the CELL sections! "//& 01370 CPSourceFileRef,& 01371 error,failure) 01372 01373 ! Handling Parallel execution 01374 CALL mp_sync(force_env%para_env%group) 01375 ! Let's transfer energy, natom, forces, virials 01376 CALL mp_sum(energies, force_env%para_env%group) 01377 CALL mp_sum(glob_natoms, force_env%para_env%group) 01378 ! Transfer forces 01379 DO iforce_eval = 1, nforce_eval 01380 ALLOCATE(global_forces(iforce_eval)%forces(3,glob_natoms(iforce_eval)),stat=stat) 01381 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01382 global_forces(iforce_eval)%forces = 0.0_dp 01383 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN 01384 IF ( force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos==& 01385 force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN 01386 ! Forces 01387 DO iparticle = 1, glob_natoms(iforce_eval) 01388 global_forces(iforce_eval)%forces(:,iparticle) = & 01389 particles(iforce_eval)%list%els(iparticle)%f 01390 END DO 01391 END IF 01392 END IF 01393 CALL mp_sum(global_forces(iforce_eval)%forces, force_env%para_env%group) 01394 !Transfer only the relevant part of the virial.. 01395 CALL mp_sum(virials(iforce_eval)%virial%pv_total, force_env%para_env%group) 01396 CALL mp_sum(virials(iforce_eval)%virial%pv_kinetic, force_env%para_env%group) 01397 CALL mp_sum(virials(iforce_eval)%virial%pv_virial, force_env%para_env%group) 01398 CALL mp_sum(virials(iforce_eval)%virial%pv_xc, force_env%para_env%group) 01399 CALL mp_sum(virials(iforce_eval)%virial%pv_fock_4c, force_env%para_env%group) 01400 CALL mp_sum(virials(iforce_eval)%virial%pv_constraint, force_env%para_env%group) 01401 !Transfer results 01402 source = 0 01403 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN 01404 IF ( force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos==& 01405 force_env%sub_force_env(iforce_eval)%force_env%para_env%source)& 01406 source=force_env%para_env%mepos 01407 ENDIF 01408 CALL mp_sum(source, force_env%para_env%group) 01409 CALL cp_results_mp_bcast(results(iforce_eval)%results, source, force_env%para_env, error) 01410 END DO 01411 01412 force_env%mixed_env%energies = energies 01413 ! Start combining the different sub_force_env 01414 CALL get_mixed_env(mixed_env=force_env%mixed_env,& 01415 mixed_energy=mixed_energy,& 01416 error=error) 01417 01418 !NB: do this for all MIXING_TYPE values, since some need it (e.g. linear mixing 01419 !NB if the first system has fewer atoms than the second) 01420 DO iparticle = 1, SIZE(particles_mix%els) 01421 particles_mix%els(iparticle)%f(:) = 0.0_dp 01422 END DO 01423 01424 CALL section_vals_val_get(mixed_section,"MIXING_TYPE",i_val=mixing_type,error=error) 01425 SELECT CASE(mixing_type) 01426 CASE(mix_linear_combination) 01427 ! Support offered only 2 force_eval 01428 CPPrecondition(nforce_eval==2,cp_failure_level,routineP,error,failure) 01429 CALL section_vals_val_get(mixed_section,"LINEAR%LAMBDA",r_val=lambda,error=error) 01430 mixed_energy%pot=lambda*energies(1) + (1.0_dp-lambda)*energies(2) 01431 ! General Mapping of forces... 01432 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01433 lambda, 1, nforce_eval, map_index, mixed_section, mapping_section, .TRUE., error) 01434 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01435 (1.0_dp-lambda), 2, nforce_eval, map_index, mixed_section, mapping_section, .FALSE., error) 01436 CASE(mix_minimum) 01437 ! Support offered only 2 force_eval 01438 CPPrecondition(nforce_eval==2,cp_failure_level,routineP,error,failure) 01439 IF (energies(1)<energies(2)) THEN 01440 mixed_energy%pot=energies(1) 01441 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01442 1.0_dp, 1, nforce_eval, map_index, mixed_section, mapping_section, .TRUE., error) 01443 ELSE 01444 mixed_energy%pot=energies(2) 01445 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01446 1.0_dp, 2, nforce_eval, map_index, mixed_section, mapping_section, .TRUE., error) 01447 ENDIF 01448 CASE(mix_coupled) 01449 ! Support offered only 2 force_eval 01450 CPPrecondition(nforce_eval==2,cp_failure_level,routineP,error,failure) 01451 CALL section_vals_val_get(mixed_section,"COUPLING%COUPLING_PARAMETER",& 01452 r_val=coupling_parameter,error=error) 01453 sd = SQRT((energies(1)-energies(2))**2+4.0_dp*coupling_parameter**2) 01454 der_1=(1.0_dp-(1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1)-energies(2)))/2.0_dp 01455 der_2=(1.0_dp+(1.0_dp/(2.0_dp*sd))*2.0_dp*(energies(1)-energies(2)))/2.0_dp 01456 mixed_energy%pot=(energies(1)+energies(2)-sd)/2.0_dp 01457 ! General Mapping of forces... 01458 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01459 der_1, 1, nforce_eval, map_index, mixed_section, mapping_section, .TRUE., error) 01460 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01461 der_2, 2, nforce_eval, map_index, mixed_section, mapping_section, .FALSE., error) 01462 CASE(mix_restrained) 01463 ! Support offered only 2 force_eval 01464 CPPrecondition(nforce_eval==2,cp_failure_level,routineP,error,failure) 01465 CALL section_vals_val_get(mixed_section,"RESTRAINT%RESTRAINT_TARGET",& 01466 r_val=restraint_target,error=error) 01467 CALL section_vals_val_get(mixed_section,"RESTRAINT%RESTRAINT_STRENGTH",& 01468 r_val=restraint_strength,error=error) 01469 mixed_energy%pot=energies(1)+restraint_strength*(energies(1)-energies(2)-restraint_target)**2 01470 der_2 = -2.0_dp*restraint_strength*(energies(1)-energies(2)-restraint_target) 01471 der_1 = 1.0_dp - der_2 01472 ! General Mapping of forces... 01473 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01474 der_1, 1, nforce_eval, map_index, mixed_section, mapping_section, .TRUE., error) 01475 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01476 der_2, 2, nforce_eval, map_index, mixed_section, mapping_section, .FALSE., error) 01477 CASE(mix_generic) 01478 ! Support any number of force_eval sections 01479 gen_section => section_vals_get_subs_vals(mixed_section,"GENERIC",error=error) 01480 CALL get_generic_info(gen_section, "MIXING_FUNCTION", coupling_function, force_env%mixed_env%par,& 01481 force_env%mixed_env%val, energies, error=error) 01482 CALL initf(1) 01483 CALL parsef(1,TRIM(coupling_function),force_env%mixed_env%par) 01484 ! Now the hardest part.. map energy with corresponding force_eval 01485 mixed_energy%pot= evalf(1,force_env%mixed_env%val) 01486 CPPrecondition(EvalErrType<=0,cp_failure_level,routineP,error,failure) 01487 CALL zero_virial(virial_mix, reset=.FALSE.) 01488 CALL cp_results_erase(results_mix, error=error) 01489 DO iforce_eval = 1, nforce_eval 01490 CALL section_vals_val_get(gen_section,"DX",r_val=dx,error=error) 01491 CALL section_vals_val_get(gen_section,"ERROR_LIMIT",r_val=lerr,error=error) 01492 dedf = evalfd(1,iforce_eval,force_env%mixed_env%val,dx,err) 01493 IF (ABS(err)>lerr) THEN 01494 WRITE(this_error,"(A,G12.6,A)")"(",err,")" 01495 WRITE(def_error,"(A,G12.6,A)")"(",lerr,")" 01496 CALL compress(this_error,.TRUE.) 01497 CALL compress(def_error,.TRUE.) 01498 CALL cp_assert(.FALSE.,cp_warning_level,-300,routineP,& 01499 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)//& 01500 ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'//& 01501 TRIM(def_error)//' .',error=error,only_ionode=.TRUE.) 01502 END IF 01503 ! General Mapping of forces... 01504 CALL mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, virials, results,& 01505 dedf, iforce_eval, nforce_eval, map_index, mixed_section, mapping_section, .FALSE., error) 01506 force_env%mixed_env%val(iforce_eval) = energies(iforce_eval) 01507 END DO 01508 ! Let's store the needed information.. 01509 force_env%mixed_env%dx = dx 01510 force_env%mixed_env%lerr= lerr 01511 force_env%mixed_env%coupling_function = TRIM(coupling_function) 01512 CALL finalizef() 01513 CASE DEFAULT 01514 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01515 END SELECT 01516 !Simply deallocate and loose the pointer references.. 01517 DO iforce_eval = 1, nforce_eval 01518 DEALLOCATE(global_forces(iforce_eval)%forces,stat=stat) 01519 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01520 CALL virial_release(virials(iforce_eval)%virial, error=error) 01521 CALL cp_result_release(results(iforce_eval)%results, error=error) 01522 END DO 01523 DEALLOCATE(global_forces, stat=stat) 01524 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01525 DEALLOCATE(subsystems, stat=stat) 01526 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01527 DEALLOCATE(particles, stat=stat) 01528 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01529 DEALLOCATE(energies, stat=stat) 01530 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01531 DEALLOCATE(glob_natoms, stat=stat) 01532 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01533 DEALLOCATE(virials, stat=stat) 01534 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01535 DEALLOCATE(results, stat=stat) 01536 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01537 ! Print Section 01538 unit_nr=cp_print_key_unit_nr(logger,mixed_section,"PRINT%DIPOLE",& 01539 extension=".data",middle_name="MIXED_DIPOLE",log_filename=.FALSE.,error=error) 01540 IF (unit_nr>0) THEN 01541 description ='[DIPOLE]' 01542 CALL get_results(results=results_mix,description=description,values=dip_mix,error=error) 01543 WRITE(unit_nr,'(/,1X,A,T48,3F11.6)')"MIXED ENV| DIPOLE ( A.U.)|",dip_mix 01544 WRITE(unit_nr,'( 1X,A,T48,3F11.6)')"MIXED ENV| DIPOLE (Debye)|",dip_mix*debye 01545 END IF 01546 CALL cp_print_key_finished_output(unit_nr,logger,mixed_section,"PRINT%DIPOLE",error=error) 01547 END SUBROUTINE mixed_energy_forces 01548 01549 END MODULE force_env_methods
1.7.3