|
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 ! ***************************************************************************** 00011 MODULE external_potential_methods 00012 USE cp_subsys_types, ONLY: cp_subsys_get,& 00013 cp_subsys_type 00014 USE f77_blas 00015 USE force_env_types, ONLY: force_env_get,& 00016 force_env_set,& 00017 force_env_type 00018 USE force_fields_util, ONLY: get_generic_info 00019 USE fparser, ONLY: evalf,& 00020 evalfd,& 00021 finalizef,& 00022 initf,& 00023 parsef 00024 USE input_section_types, ONLY: section_vals_get,& 00025 section_vals_get_subs_vals,& 00026 section_vals_type,& 00027 section_vals_val_get 00028 USE kinds, ONLY: default_path_length,& 00029 default_string_length,& 00030 dp 00031 USE memory_utilities, ONLY: reallocate 00032 USE particle_list_types, ONLY: particle_list_type 00033 USE string_utilities, ONLY: compress 00034 USE timings, ONLY: timeset,& 00035 timestop 00036 #include "cp_common_uses.h" 00037 00038 IMPLICIT NONE 00039 00040 PRIVATE 00041 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'external_potential_methods' 00042 PUBLIC :: add_external_potential 00043 00044 CONTAINS 00045 00046 ! ***************************************************************************** 00051 SUBROUTINE add_external_potential(force_env, error) 00052 TYPE(force_env_type), POINTER :: force_env 00053 TYPE(cp_error_type), INTENT(inout) :: error 00054 00055 CHARACTER(len=*), PARAMETER :: routineN = 'add_external_potential', 00056 routineP = moduleN//':'//routineN 00057 00058 CHARACTER(LEN=default_path_length) :: coupling_function 00059 CHARACTER(LEN=default_string_length) :: def_error, this_error 00060 CHARACTER(LEN=default_string_length), 00061 DIMENSION(:), POINTER :: my_par 00062 INTEGER :: a_var, handle, i, iatom, j, 00063 k, n_var, natom, rep, stat 00064 INTEGER, DIMENSION(:), POINTER :: iatms, nparticle 00065 LOGICAL :: failure, useall 00066 REAL(KIND=dp) :: dedf, dx, energy, err, lerr 00067 REAL(KIND=dp), DIMENSION(:), POINTER :: my_val 00068 TYPE(cp_logger_type), POINTER :: logger 00069 TYPE(cp_subsys_type), POINTER :: subsys 00070 TYPE(particle_list_type), POINTER :: particles 00071 TYPE(section_vals_type), POINTER :: ext_pot_section 00072 00073 failure = .FALSE. 00074 useall = .FALSE. 00075 CALL timeset(routineN,handle) 00076 NULLIFY(my_par, my_val, logger, subsys, particles, ext_pot_section, nparticle) 00077 ext_pot_section => section_vals_get_subs_vals(force_env%force_env_section,& 00078 "EXTERNAL_POTENTIAL",error=error) 00079 CALL section_vals_get(ext_pot_section,n_repetition=n_var,error=error) 00080 DO rep=1, n_var 00081 natom = 0 00082 logger => cp_error_get_logger(error) 00083 CALL section_vals_val_get(ext_pot_section,"DX",r_val=dx,i_rep_section=rep,error=error) 00084 CALL section_vals_val_get(ext_pot_section,"ERROR_LIMIT",r_val=lerr,i_rep_section=rep,error=error) 00085 CALL get_generic_info(ext_pot_section, "FUNCTION", coupling_function, my_par, my_val,& 00086 input_variables=(/"X","Y","Z"/), i_rep_sec=rep,error=error) 00087 CALL initf(1) 00088 CALL parsef(1,TRIM(coupling_function),my_par) 00089 00090 ! Apply potential on all atoms, computing energy and forces 00091 NULLIFY(particles, subsys) 00092 CALL force_env_get(force_env, subsys=subsys, error=error) 00093 CALL cp_subsys_get(subsys, particles=particles, error=error) 00094 CALL force_env_get(force_env, additional_potential=energy, error=error) 00095 CALL section_vals_val_get(ext_pot_section,"ATOMS_LIST",n_rep_val=a_var,i_rep_section=rep,error=error) 00096 DO k = 1, a_var 00097 CALL section_vals_val_get(ext_pot_section,"ATOMS_LIST",i_rep_val=k,i_vals=iatms,i_rep_section=rep,& 00098 error=error) 00099 CALL reallocate(nparticle,1, natom+SIZE(iatms)) 00100 nparticle(natom+1:natom+SIZE(iatms)) = iatms 00101 natom = natom + SIZE(iatms) 00102 END DO 00103 IF (a_var==0) THEN 00104 natom = particles%n_els 00105 useall = .TRUE. 00106 END IF 00107 DO i = 1, natom 00108 IF (useall) THEN 00109 iatom = i 00110 ELSE 00111 iatom = nparticle(i) 00112 END IF 00113 my_val(1)=particles%els(iatom)%r(1) 00114 my_val(2)=particles%els(iatom)%r(2) 00115 my_val(3)=particles%els(iatom)%r(3) 00116 00117 energy = energy + evalf(1,my_val) 00118 DO j = 1, 3 00119 dedf = evalfd(1,j,my_val,dx,err) 00120 IF (ABS(err)>lerr) THEN 00121 WRITE(this_error,"(A,G12.6,A)")"(",err,")" 00122 WRITE(def_error,"(A,G12.6,A)")"(",lerr,")" 00123 CALL compress(this_error,.TRUE.) 00124 CALL compress(def_error,.TRUE.) 00125 CALL cp_assert(.FALSE.,cp_warning_level,-300,routineP,& 00126 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)//& 00127 ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'//& 00128 TRIM(def_error)//' .',error=error,only_ionode=.TRUE.) 00129 END IF 00130 particles%els(iatom)%f(j)=particles%els(iatom)%f(j)-dedf 00131 END DO 00132 END DO 00133 CALL force_env_set(force_env, additional_potential=energy, error=error) 00134 DEALLOCATE(my_par,stat=stat) 00135 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00136 DEALLOCATE(my_val,stat=stat) 00137 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00138 IF (a_var/=0) THEN 00139 DEALLOCATE (nparticle,STAT=stat) 00140 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00141 END IF 00142 CALL finalizef() 00143 END DO 00144 CALL timestop(handle) 00145 END SUBROUTINE add_external_potential 00146 00147 END MODULE external_potential_methods
1.7.3