|
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 MODULE optimize_input 00006 USE cell_types, ONLY: parse_cell_line 00007 USE cp_external_control, ONLY: external_control 00008 USE cp_output_handling, ONLY: cp_add_iter_level,& 00009 cp_iterate,& 00010 cp_print_key_finished_output,& 00011 cp_print_key_unit_nr,& 00012 cp_rm_iter_level 00013 USE cp_para_types, ONLY: cp_para_env_type 00014 USE cp_parser_methods, ONLY: parser_read_line 00015 USE cp_parser_types, ONLY: cp_parser_type,& 00016 parser_create,& 00017 parser_release 00018 USE f77_blas 00019 USE f77_interface, ONLY: calc_force,& 00020 create_force_env,& 00021 destroy_force_env,& 00022 set_cell 00023 USE input_constants, ONLY: opt_force_matching 00024 USE input_cp2k_restarts, ONLY: write_restart_header 00025 USE input_section_types, ONLY: section_vals_get,& 00026 section_vals_get_subs_vals,& 00027 section_vals_type,& 00028 section_vals_val_get,& 00029 section_vals_val_set,& 00030 section_vals_write 00031 USE kinds, ONLY: default_path_length,& 00032 default_string_length,& 00033 dp 00034 USE machine, ONLY: m_flush,& 00035 m_walltime 00036 USE memory_utilities, ONLY: reallocate 00037 USE message_passing, ONLY: mp_bcast,& 00038 mp_comm_free,& 00039 mp_comm_split,& 00040 mp_comm_split_direct,& 00041 mp_environ,& 00042 mp_sum 00043 USE parallel_rng_types, ONLY: UNIFORM,& 00044 create_rng_stream,& 00045 delete_rng_stream,& 00046 next_random_number,& 00047 rng_stream_type 00048 USE physcon, ONLY: bohr 00049 USE powell, ONLY: opt_state_type,& 00050 powell_optimize 00051 USE timings, ONLY: timeset,& 00052 timestop 00053 #include "cp_common_uses.h" 00054 00055 IMPLICIT NONE 00056 PRIVATE 00057 00058 PUBLIC :: run_optimize_input 00059 00060 TYPE fm_env_type 00061 CHARACTER(LEN=default_path_length) :: optimize_file_name 00062 00063 CHARACTER(LEN=default_path_length) :: ref_traj_file_name 00064 CHARACTER(LEN=default_path_length) :: ref_force_file_name 00065 CHARACTER(LEN=default_path_length) :: ref_cell_file_name 00066 00067 INTEGER :: group_size 00068 00069 REAL(KIND=dp) :: energy_weight 00070 REAL(KIND=dp) :: shift_mm 00071 REAL(KIND=dp) :: shift_qm 00072 LOGICAL :: shift_average 00073 INTEGER :: frame_start,frame_stop,frame_stride,frame_count 00074 END TYPE 00075 00076 TYPE variable_type 00077 CHARACTER(LEN=default_string_length) :: label 00078 REAL(KIND=dp) :: value 00079 LOGICAL :: fixed 00080 END TYPE 00081 00082 TYPE oi_env_type 00083 INTEGER :: method 00084 INTEGER :: seed 00085 CHARACTER(LEN=default_path_length) :: project_name 00086 TYPE(fm_env_type) :: fm_env 00087 TYPE(variable_type), DIMENSION(:), ALLOCATABLE :: variables 00088 REAL(KIND=dp) :: rhobeg,rhoend 00089 INTEGER :: maxfun 00090 INTEGER :: iter_start_val 00091 REAL(KIND=dp) :: randomize_variables 00092 REAL(KIND=dp) :: start_time,target_time 00093 END TYPE 00094 00095 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_input' 00096 00097 CONTAINS 00098 00099 ! ***************************************************************************** 00103 SUBROUTINE run_optimize_input(root_section,para_env, error) 00104 TYPE(section_vals_type), POINTER :: root_section 00105 TYPE(cp_para_env_type), POINTER :: para_env 00106 TYPE(cp_error_type), INTENT(INOUT) :: error 00107 00108 CHARACTER(len=*), PARAMETER :: routineN = 'run_optimize_input', 00109 routineP = moduleN//':'//routineN 00110 00111 INTEGER :: handle, i_var 00112 LOGICAL :: failure 00113 REAL(KIND=dp) :: random_number, seed(3,2) 00114 TYPE(oi_env_type) :: oi_env 00115 TYPE(rng_stream_type), POINTER :: rng_stream 00116 00117 failure=.FALSE. 00118 00119 CALL timeset(routineN,handle) 00120 00121 oi_env%start_time=m_walltime() 00122 00123 CALL parse_input(oi_env,root_section,error) 00124 00125 ! if we have been asked to randomize the variables, we do this. 00126 IF (oi_env%randomize_variables.NE.0.0_dp) THEN 00127 NULLIFY(rng_stream) 00128 seed=REAL(oi_env%seed,KIND=dp) 00129 CALL create_rng_stream(rng_stream,"run_optimize_input",distribution_type=UNIFORM,seed=seed,error=error) 00130 DO i_var=1,SIZE(oi_env%variables,1) 00131 IF (.NOT.oi_env%variables(i_var)%fixed) THEN 00132 ! change with a random percentage the variable 00133 random_number = next_random_number(rng_stream,error=error) 00134 oi_env%variables(i_var)%value=oi_env%variables(i_var)%value* & 00135 (1.0_dp+(2*random_number-1.0_dp)*oi_env%randomize_variables/100.0_dp) 00136 ENDIF 00137 ENDDO 00138 CALL delete_rng_stream(rng_stream,error=error) 00139 ENDIF 00140 00141 ! proceed to actual methods 00142 SELECT CASE(oi_env%method) 00143 CASE(opt_force_matching) 00144 CALL force_matching(oi_env,root_section,para_env,error) 00145 CASE DEFAULT 00146 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00147 END SELECT 00148 00149 CALL timestop(handle) 00150 00151 END SUBROUTINE run_optimize_input 00152 00153 ! ***************************************************************************** 00157 SUBROUTINE force_matching(oi_env,root_section,para_env, error) 00158 TYPE(oi_env_type) :: oi_env 00159 TYPE(section_vals_type), POINTER :: root_section 00160 TYPE(cp_para_env_type), POINTER :: para_env 00161 TYPE(cp_error_type), INTENT(INOUT) :: error 00162 00163 CHARACTER(len=*), PARAMETER :: routineN = 'force_matching', 00164 routineP = moduleN//':'//routineN 00165 00166 CHARACTER(len=default_path_length) :: input_path, output_path 00167 CHARACTER(len=default_string_length), 00168 ALLOCATABLE, DIMENSION(:, :) :: initial_variables 00169 INTEGER :: color, energies_unit, handle, history_unit, i_atom, i_el, 00170 i_frame, i_free_var, i_var, ierr, mepos_master, mepos_slave, 00171 mpi_comm_master, mpi_comm_slave, mpi_comm_slave_primus, n_atom, n_el, 00172 n_frames, n_free_var, n_groups, n_var, new_env_id, num_pe_master, 00173 num_pe_slave, output_unit, restart_unit, state 00174 INTEGER, ALLOCATABLE, DIMENSION(:) :: free_var_index 00175 INTEGER, DIMENSION(:), POINTER :: group_distribution 00176 LOGICAL :: should_stop 00177 REAL(KIND=dp) :: e1, e2, e3, e4, e_pot, 00178 energy_weight, re, rf, 00179 shift_mm, shift_qm, t1, t2, 00180 t3, t4, t5 00181 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: force, free_vars, pos 00182 REAL(KIND=dp), DIMENSION(:), POINTER :: energy_traj, 00183 energy_traj_read, energy_var 00184 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cell_traj, cell_traj_read, 00185 force_traj, force_traj_read, force_var, pos_traj, pos_traj_read 00186 TYPE(cp_logger_type), POINTER :: logger 00187 TYPE(opt_state_type) :: ostate 00188 TYPE(section_vals_type), POINTER :: oi_section, variable_section 00189 00190 CALL timeset(routineN,handle) 00191 00192 logger=>cp_error_get_logger(error) 00193 CALL cp_add_iter_level(logger%iter_info,"POWELL_OPT",error=error) 00194 output_unit = cp_logger_get_default_io_unit(logger) 00195 00196 IF (output_unit>0) THEN 00197 WRITE(output_unit,'(T2,A)') 'FORCE_MATCHING| good morning....' 00198 ENDIF 00199 00200 ! do IO of ref traj / frc / cell 00201 NULLIFY(cell_traj_read,force_traj_read,pos_traj_read,energy_traj_read) 00202 CALL read_reference_data(oi_env,para_env,force_traj_read,pos_traj_read,energy_traj_read,cell_traj_read,error) 00203 n_atom=SIZE(pos_traj_read,2) 00204 00205 ! adjust read data with respect to start/stop/stride 00206 IF (oi_env%fm_env%frame_stop<0) oi_env%fm_env%frame_stop=SIZE(pos_traj_read,3) 00207 00208 IF (oi_env%fm_env%frame_count>0) THEN 00209 oi_env%fm_env%frame_stride=(oi_env%fm_env%frame_stop-oi_env%fm_env%frame_start+1+ & 00210 oi_env%fm_env%frame_count-1)/(oi_env%fm_env%frame_count) 00211 ENDIF 00212 n_frames=(oi_env%fm_env%frame_stop-oi_env%fm_env%frame_start+oi_env%fm_env%frame_stride)/oi_env%fm_env%frame_stride 00213 00214 ALLOCATE(force_traj(3,n_atom,n_frames),pos_traj(3,n_atom,n_frames),energy_traj(n_frames)) 00215 IF (ASSOCIATED(cell_traj_read)) ALLOCATE(cell_traj(3,3,n_frames)) 00216 00217 n_frames=0 00218 DO i_frame=oi_env%fm_env%frame_start,oi_env%fm_env%frame_stop,oi_env%fm_env%frame_stride 00219 n_frames=n_frames+1 00220 force_traj(:,:,n_frames)=force_traj_read(:,:,i_frame) 00221 pos_traj(:,:,n_frames)=pos_traj_read(:,:,i_frame) 00222 energy_traj(n_frames)=energy_traj_read(i_frame) 00223 IF (ASSOCIATED(cell_traj)) cell_traj(:,:,n_frames)=cell_traj_read(:,:,i_frame) 00224 ENDDO 00225 DEALLOCATE(force_traj_read,pos_traj_read,energy_traj_read) 00226 IF (ASSOCIATED(cell_traj_read)) DEALLOCATE(cell_traj_read) 00227 00228 n_el=3*n_atom 00229 ALLOCATE(pos(n_el), force(n_el)) 00230 ALLOCATE(energy_var(n_frames),force_var(3,n_atom,n_frames)) 00231 00232 00233 ! split the para_env in a set of sub_para_envs that will do the force_env communications 00234 mpi_comm_master=para_env%group 00235 num_pe_master =para_env%num_pe 00236 mepos_master =para_env%mepos 00237 ALLOCATE(group_distribution(0:num_pe_master-1)) 00238 IF (oi_env%fm_env%group_size>para_env%num_pe) oi_env%fm_env%group_size=para_env%num_pe 00239 00240 CALL mp_comm_split(mpi_comm_master,mpi_comm_slave,n_groups,group_distribution,subgroup_min_size=oi_env%fm_env%group_size) 00241 CALL mp_environ(num_pe_slave,mepos_slave,mpi_comm_slave) 00242 color=0 00243 IF (mepos_slave==0) color=1 00244 CALL mp_comm_split_direct(mpi_comm_master,mpi_comm_slave_primus,color) 00245 00246 ! assign initial variables 00247 n_var=SIZE(oi_env%variables,1) 00248 ALLOCATE(initial_variables(2,n_var)) 00249 n_free_var=0 00250 DO i_var=1,n_var 00251 initial_variables(1,i_var)=oi_env%variables(i_var)%label 00252 WRITE(initial_variables(2,i_var),*) oi_env%variables(i_var)%value 00253 IF (.NOT.oi_env%variables(i_var)%fixed) n_free_var=n_free_var+1 00254 ENDDO 00255 ALLOCATE(free_vars(n_free_var),free_var_index(n_free_var)) 00256 i_free_var=0 00257 DO i_var=1,n_var 00258 IF (.NOT.oi_env%variables(i_var)%fixed) THEN 00259 i_free_var=i_free_var+1 00260 free_var_index(i_free_var)=i_var 00261 free_vars(i_free_var)=oi_env%variables(free_var_index(i_free_var))%value 00262 ENDIF 00263 ENDDO 00264 00265 ! create input and output file names. 00266 input_path=oi_env%fm_env%optimize_file_name 00267 WRITE(output_path,'(A,I0,A)') TRIM(oi_env%project_name)//"-worker-",group_distribution(mepos_master),".out" 00268 00269 ! initialize the powell optimizer 00270 energy_weight=oi_env%fm_env%energy_weight 00271 shift_mm=oi_env%fm_env%shift_mm 00272 shift_qm=oi_env%fm_env%shift_qm 00273 00274 IF (para_env%mepos==para_env%source) THEN 00275 ostate%nf = 0 00276 ostate%nvar = n_free_var 00277 ostate%rhoend = oi_env%rhoend 00278 ostate%rhobeg = oi_env%rhobeg 00279 ostate%maxfun = oi_env%maxfun 00280 ostate%iprint = 1 00281 ostate%unit = output_unit 00282 ostate%state = 0 00283 ENDIF 00284 00285 IF (output_unit>0) THEN 00286 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| number of atoms per frame ',n_atom 00287 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| number of frames ',n_frames 00288 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| number of parallel groups ',n_groups 00289 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| number of variables ',n_var 00290 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| number of free variables ',n_free_var 00291 WRITE(output_unit,'(T2,A,A)') 'FORCE_MATCHING| optimize file name ',TRIM(input_path) 00292 WRITE(output_unit,'(T2,A,T60,F20.12)') 'FORCE_MATCHING| accuracy',ostate%rhoend 00293 WRITE(output_unit,'(T2,A,T60,F20.12)') 'FORCE_MATCHING| step size',ostate%rhobeg 00294 WRITE(output_unit,'(T2,A,T60,I20)') 'FORCE_MATCHING| max function evaluation',ostate%maxfun 00295 WRITE(output_unit,'(T2,A,T60,L20)') 'FORCE_MATCHING| shift average',oi_env%fm_env%shift_average 00296 WRITE(output_unit,'(T2,A)') 'FORCE_MATCHING| initial values:' 00297 DO i_var=1,n_var 00298 WRITE(output_unit,'(T2,A,1X,E28.16)') TRIM(oi_env%variables(i_var)%label),oi_env%variables(i_var)%value 00299 ENDDO 00300 WRITE(output_unit,'(T2,A)') 'FORCE_MATCHING| switching to POWELL optimization of the free parameters' 00301 WRITE(output_unit,'()') 00302 WRITE(output_unit,'(T2,A20,A20,A11,A11)') 'iteration number','function value','time','time Force' 00303 CALL m_flush(output_unit) 00304 ENDIF 00305 00306 00307 t1 = m_walltime() 00308 00309 DO 00310 00311 00312 ! globalize the state 00313 IF (para_env%mepos==para_env%source) state=ostate%state 00314 CALL mp_bcast(state,para_env%source,para_env%group) 00315 00316 ! if required get the energy of this set of params 00317 IF (state == 2) THEN 00318 00319 CALL cp_iterate(logger%iter_info,last=.FALSE.,error=error) 00320 00321 ! create a new force env, updating the free vars as needed 00322 DO i_free_var=1,n_free_var 00323 WRITE(initial_variables(2,free_var_index(i_free_var)),*) free_vars(i_free_var) 00324 oi_env%variables(free_var_index(i_free_var))%value=free_vars(i_free_var) 00325 ENDDO 00326 00327 ierr=0 00328 CALL create_force_env(new_env_id=new_env_id,input_path=input_path, output_path=output_path, & 00329 mpi_comm=mpi_comm_slave,initial_variables=initial_variables,ierr=ierr) 00330 00331 ! set to zero initialy, for easier mp_summing 00332 energy_var=0.0_dp 00333 force_var=0.0_dp 00334 00335 ! compute energies and forces for all frames, doing the work on a slave sub group based on round robin 00336 t5 = 0.0_dp 00337 DO i_frame=group_distribution(mepos_master)+1,n_frames,n_groups 00338 00339 ! set new cell if needed 00340 IF (ASSOCIATED(cell_traj)) THEN 00341 CALL set_cell(env_id=new_env_id, new_cell=cell_traj(:,:,i_frame), ierr=ierr) 00342 ENDIF 00343 00344 ! copy pos from ref 00345 i_el=0 00346 DO i_atom=1,n_atom 00347 pos(i_el+1)=pos_traj(1,i_atom,i_frame) 00348 pos(i_el+2)=pos_traj(2,i_atom,i_frame) 00349 pos(i_el+3)=pos_traj(3,i_atom,i_frame) 00350 i_el=i_el+3 00351 ENDDO 00352 00353 ! evaluate energy/force with new pos 00354 t3 = m_walltime() 00355 CALL calc_force(env_id=new_env_id,pos=pos,n_el_pos=n_el,e_pot=e_pot,force=force,n_el_force=n_el,ierr=ierr) 00356 t4 = m_walltime() 00357 t5 = t5 + t4-t3 00358 00359 ! copy force and energy in place 00360 energy_var(i_frame)=e_pot 00361 i_el=0 00362 DO i_atom=1,n_atom 00363 force_var(1,i_atom,i_frame)=force(i_el+1) 00364 force_var(2,i_atom,i_frame)=force(i_el+2) 00365 force_var(3,i_atom,i_frame)=force(i_el+3) 00366 i_el=i_el+3 00367 ENDDO 00368 00369 ENDDO 00370 00371 ! clean up force env, get ready for the next round 00372 CALL destroy_force_env(env_id=new_env_id,ierr=ierr) 00373 00374 ! get data everywhere on the master group, we could reduce the amount of data by reducing to partial RMSD first 00375 ! furthermore, we should only do this operation among the masters of the slave group 00376 IF (mepos_slave==0) THEN 00377 CALL mp_sum(energy_var,mpi_comm_slave_primus) 00378 CALL mp_sum(force_var,mpi_comm_slave_primus) 00379 ENDIF 00380 00381 ! now evaluate the target function to be minimized (only valid on mepos_slave==0) 00382 IF (para_env%mepos==para_env%source) THEN 00383 rf = SQRT(SUM((force_var-force_traj)**2)/(REAL(n_frames,KIND=dp) * REAL(n_atom,KIND=dp))) 00384 IF (oi_env%fm_env%shift_average) THEN 00385 shift_mm=SUM(energy_var)/n_frames 00386 shift_qm=SUM(energy_traj)/n_frames 00387 ENDIF 00388 re = SQRT(SUM(((energy_var-shift_mm)-(energy_traj-shift_qm))**2)/n_frames) 00389 ostate%f= energy_weight * re + rf 00390 t2 = m_walltime() 00391 WRITE(output_unit,'(T2,I20,F20.12,F11.3,F11.3)') oi_env%iter_start_val + ostate%nf, ostate%f, t2-t1, t5 00392 CALL m_flush(output_unit) 00393 t1 = m_walltime() 00394 ENDIF 00395 00396 ! the history file with the trajectory of the parameters 00397 history_unit=cp_print_key_unit_nr(logger,root_section,"OPTIMIZE_INPUT%HISTORY", & 00398 extension=".dat",error=error) 00399 IF (history_unit>0) THEN 00400 WRITE (UNIT=history_unit,FMT="(I20,F20.12,1000F20.12)") oi_env%iter_start_val + ostate%nf, ostate%f, free_vars 00401 END IF 00402 CALL cp_print_key_finished_output(history_unit,logger,root_section,"OPTIMIZE_INPUT%HISTORY", & 00403 error=error) 00404 00405 ! the energy profile for all frames 00406 energies_unit=cp_print_key_unit_nr(logger,root_section,"OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_ENERGIES", & 00407 file_position="REWIND",extension=".dat",error=error) 00408 IF (energies_unit>0) THEN 00409 WRITE (UNIT=energies_unit,FMT="(A20,A20,A20,A20)") "#frame","ref","fit","diff" 00410 DO i_frame=1,n_frames 00411 e1=energy_traj(i_frame)-shift_qm 00412 e2=energy_var(i_frame)-shift_mm 00413 WRITE (UNIT=energies_unit,FMT="(I20,F20.12,F20.12,F20.12)") i_frame,e1,e2,e1-e2 00414 ENDDO 00415 END IF 00416 CALL cp_print_key_finished_output(energies_unit,logger,root_section,"OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_ENERGIES", & 00417 error=error) 00418 00419 ! the force profile for all frames 00420 energies_unit=cp_print_key_unit_nr(logger,root_section,"OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_FORCES", & 00421 file_position="REWIND",extension=".dat",error=error) 00422 IF (energies_unit>0) THEN 00423 WRITE (UNIT=energies_unit,FMT="(A20,A20,A20,A20)") "#frame","normalized diff","diff","ref","ref sum" 00424 DO i_frame=1,n_frames 00425 e1=SQRT(SUM((force_var(:,:,i_frame)-force_traj(:,:,i_frame))**2)) 00426 e2=SQRT(SUM((force_traj(:,:,i_frame))**2)) 00427 e3=SQRT(SUM(SUM(force_traj(:,:,i_frame),DIM=2)**2)) 00428 e4=SQRT(SUM(SUM(force_var(:,:,i_frame),DIM=2)**2)) 00429 WRITE (UNIT=energies_unit,FMT="(I20,F20.12,F20.12,F20.12,2F20.12)") i_frame,e1/e2,e1,e2,e3,e4 00430 ENDDO 00431 END IF 00432 CALL cp_print_key_finished_output(energies_unit,logger,root_section,"OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_FORCES", & 00433 error=error) 00434 00435 00436 ! a restart file with the current values of the parameters 00437 restart_unit = cp_print_key_unit_nr(logger,root_section,"OPTIMIZE_INPUT%RESTART", extension=".restart",& 00438 file_position="REWIND", do_backup=.TRUE., error=error) 00439 IF (restart_unit>0) THEN 00440 oi_section => section_vals_get_subs_vals(root_section,"OPTIMIZE_INPUT",error=error) 00441 CALL section_vals_val_set(oi_section,"ITER_START_VAL",i_val=oi_env%iter_start_val + ostate%nf,error=error) 00442 variable_section => section_vals_get_subs_vals(oi_section,"VARIABLE",error=error) 00443 DO i_free_var=1,n_free_var 00444 CALL section_vals_val_set(variable_section,"VALUE",i_rep_section=free_var_index(i_free_var), & 00445 r_val=free_vars(i_free_var),error=error) 00446 ENDDO 00447 CALL write_restart_header(restart_unit, error) 00448 CALL section_vals_write(root_section,unit_nr=restart_unit,hide_root=.TRUE., error=error) 00449 ENDIF 00450 CALL cp_print_key_finished_output(restart_unit,logger,root_section,"OPTIMIZE_INPUT%RESTART", error=error) 00451 00452 ENDIF 00453 00454 IF ( state == -1 ) EXIT 00455 00456 CALL external_control(should_stop,"OPTIMIZE_INPUT",target_time=oi_env%target_time,start_time=oi_env%start_time,error=error) 00457 00458 IF (should_stop) EXIT 00459 00460 ! do a powell step if needed 00461 IF (para_env%mepos==para_env%source) THEN 00462 CALL powell_optimize (ostate%nvar, free_vars , ostate) 00463 ENDIF 00464 CALL mp_bcast(free_vars,para_env%source,para_env%group) 00465 00466 ENDDO 00467 00468 ! finally, get the best set of variables 00469 IF (para_env%mepos==para_env%source) THEN 00470 ostate%state = 8 00471 CALL powell_optimize (ostate%nvar, free_vars , ostate) 00472 ENDIF 00473 CALL mp_bcast(free_vars,para_env%source,para_env%group) 00474 DO i_free_var=1,n_free_var 00475 WRITE(initial_variables(2,free_var_index(i_free_var)),*) free_vars(i_free_var) 00476 oi_env%variables(free_var_index(i_free_var))%value=free_vars(i_free_var) 00477 ENDDO 00478 IF (para_env%mepos==para_env%source) THEN 00479 WRITE(output_unit,'(T2,A)') '' 00480 WRITE(output_unit,'(T2,A,T60,F20.12)') 'FORCE_MATCHING| optimal function value found so far:',ostate%fopt 00481 WRITE(output_unit,'(T2,A)') 'FORCE_MATCHING| optimal variables found so far:' 00482 DO i_var=1,n_var 00483 WRITE(output_unit,'(T2,A,1X,E28.16)') TRIM(oi_env%variables(i_var)%label),oi_env%variables(i_var)%value 00484 ENDDO 00485 ENDIF 00486 00487 CALL cp_rm_iter_level(logger%iter_info,"POWELL_OPT",error=error) 00488 00489 ! deallocate for cleanup 00490 IF (ASSOCIATED(cell_traj)) DEALLOCATE(cell_traj) 00491 DEALLOCATE(pos,force,force_traj, pos_traj, force_var) 00492 DEALLOCATE(group_distribution,energy_traj, energy_var) 00493 CALL mp_comm_free(mpi_comm_slave) 00494 CALL mp_comm_free(mpi_comm_slave_primus) 00495 00496 CALL timestop(handle) 00497 00498 END SUBROUTINE force_matching 00499 00500 ! ***************************************************************************** 00508 SUBROUTINE read_reference_data(oi_env,para_env,force_traj,pos_traj,energy_traj,cell_traj,error) 00509 TYPE(oi_env_type) :: oi_env 00510 TYPE(cp_para_env_type), POINTER :: para_env 00511 REAL(KIND=dp), DIMENSION(:, :, :), 00512 POINTER :: force_traj, pos_traj 00513 REAL(KIND=dp), DIMENSION(:), POINTER :: energy_traj 00514 REAL(KIND=dp), DIMENSION(:, :, :), 00515 POINTER :: cell_traj 00516 TYPE(cp_error_type), INTENT(INOUT) :: error 00517 00518 CHARACTER(len=*), PARAMETER :: routineN = 'read_reference_data', 00519 routineP = moduleN//':'//routineN 00520 00521 CHARACTER(len=default_path_length) :: filename, title 00522 CHARACTER(len=default_string_length) :: AA 00523 INTEGER :: cell_itimes, handle, i, 00524 iframe, n_frames, 00525 n_frames_current, nread, 00526 trj_itimes 00527 LOGICAL :: at_end, test_ok 00528 REAL(KIND=dp) :: cell_time, trj_epot, 00529 trj_time, vec(3), vol 00530 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: force 00531 TYPE(cp_parser_type), POINTER :: local_parser 00532 00533 CALL timeset(routineN,handle) 00534 00535 ! do IO of ref traj / frc / cell 00536 00537 ! trajectory 00538 n_frames=0 00539 n_frames_current=0 00540 NULLIFY(local_parser,pos_traj,energy_traj,force_traj) 00541 filename=oi_env%fm_env%ref_traj_file_name 00542 CALL cp_assert(filename.NE."",cp_fatal_level,cp_assertion_failed,routineP,& 00543 "The reference trajectory file name is empty"//& 00544 CPSourceFileRef,& 00545 only_ionode=.TRUE.) 00546 CALL parser_create(local_parser,filename,para_env=para_env,error=error) 00547 DO 00548 CALL parser_read_line(local_parser,1,at_end=at_end,error=error) 00549 IF (at_end) EXIT 00550 READ(local_parser%input_line,FMT="(I8)") nread 00551 n_frames=n_frames+1 00552 00553 IF (n_frames>n_frames_current) THEN 00554 n_frames_current=5*(n_frames_current+10)/3 00555 CALL reallocate(pos_traj,1,3,1,nread,1,n_frames_current) 00556 ENDIF 00557 00558 ! title line 00559 CALL parser_read_line(local_parser,1,error=error) 00560 00561 ! actual coordinates 00562 DO i = 1,nread 00563 CALL parser_read_line(local_parser,1,error=error) 00564 READ(local_parser%input_line(1:LEN_TRIM(local_parser%input_line)),*) AA,vec 00565 pos_traj(:,i,n_frames)=vec*bohr 00566 END DO 00567 00568 ENDDO 00569 CALL parser_release(local_parser,error=error) 00570 00571 n_frames_current=n_frames 00572 CALL reallocate(energy_traj,1,n_frames_current) 00573 CALL reallocate(force_traj,1,3,1,nread,1,n_frames_current) 00574 CALL reallocate(pos_traj,1,3,1,nread,1,n_frames_current) 00575 00576 ! now force reference trajectory 00577 filename=oi_env%fm_env%ref_force_file_name 00578 CALL cp_assert(filename.NE."",cp_fatal_level,cp_assertion_failed,routineP,& 00579 "The reference force file name is empty"//& 00580 CPSourceFileRef,& 00581 only_ionode=.TRUE.) 00582 CALL parser_create(local_parser,filename,para_env=para_env,error=error) 00583 DO iframe=1,n_frames 00584 CALL parser_read_line(local_parser,1,error=error) 00585 READ(local_parser%input_line,FMT="(I8)") nread 00586 00587 ! title line 00588 test_ok = .FALSE. 00589 CALL parser_read_line(local_parser,1,error=error) 00590 READ(local_parser%input_line,FMT="(T6,I8,T23,F12.3,T41,F20.10)",ERR=999) trj_itimes, trj_time, trj_epot 00591 test_ok = .TRUE. 00592 999 CONTINUE 00593 IF (.NOT. test_ok) THEN 00594 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00595 "Could not parse the title line of the trajectory file"//& 00596 CPSourceFileRef,& 00597 only_ionode=.TRUE.) 00598 END IF 00599 energy_traj(iframe)=trj_epot 00600 00601 ! actual forces, in a.u. 00602 DO i = 1,nread 00603 CALL parser_read_line(local_parser,1,error=error) 00604 READ(local_parser%input_line(1:LEN_TRIM(local_parser%input_line)),*) AA,vec 00605 force_traj(:,i,iframe)=vec 00606 END DO 00607 ENDDO 00608 CALL parser_release(local_parser,error=error) 00609 00610 00611 ! and cell, which is optional 00612 NULLIFY(cell_traj) 00613 filename=oi_env%fm_env%ref_cell_file_name 00614 IF (filename.NE."") THEN 00615 CALL parser_create(local_parser,filename,para_env=para_env,error=error) 00616 ALLOCATE(cell_traj(3,3,n_frames)) 00617 DO iframe=1,n_frames 00618 CALL parser_read_line(local_parser,1,error=error) 00619 CALL parse_cell_line(local_parser%input_line, cell_itimes, cell_time, cell_traj(:,:,iframe), vol, error) 00620 ENDDO 00621 CALL parser_release(local_parser,error=error) 00622 ENDIF 00623 00624 CALL timestop(handle) 00625 00626 END SUBROUTINE read_reference_data 00627 00628 ! ***************************************************************************** 00633 SUBROUTINE parse_input(oi_env,root_section,error) 00634 TYPE(oi_env_type) :: oi_env 00635 TYPE(section_vals_type), POINTER :: root_section 00636 TYPE(cp_error_type), INTENT(INOUT) :: error 00637 00638 CHARACTER(len=*), PARAMETER :: routineN = 'parse_input', 00639 routineP = moduleN//':'//routineN 00640 00641 INTEGER :: handle, ivar, method, n_var 00642 LOGICAL :: explicit, failure 00643 TYPE(section_vals_type), POINTER :: fm_section, oi_section, 00644 variable_section 00645 00646 CALL timeset(routineN,handle) 00647 failure=.FALSE. 00648 00649 CALL section_vals_val_get(root_section,"GLOBAL%PROJECT",c_val=oi_env%project_name,error=error) 00650 CALL section_vals_val_get(root_section,"GLOBAL%SEED",i_val=oi_env%seed,error=error) 00651 CALL section_vals_val_get(root_section,"GLOBAL%WALLTIME",r_val=oi_env%target_time,error=error) 00652 00653 oi_section => section_vals_get_subs_vals(root_section,"OPTIMIZE_INPUT",error=error) 00654 variable_section => section_vals_get_subs_vals(oi_section,"VARIABLE",error=error) 00655 00656 CALL section_vals_val_get(oi_section,"ACCURACY", r_val=oi_env%rhoend, error=error) 00657 CALL section_vals_val_get(oi_section,"STEP_SIZE", r_val=oi_env%rhobeg, error=error) 00658 CALL section_vals_val_get(oi_section,"MAX_FUN", i_val=oi_env%maxfun, error=error) 00659 CALL section_vals_val_get(oi_section,"ITER_START_VAL", i_val=oi_env%iter_start_val, error=error) 00660 CALL section_vals_val_get(oi_section,"RANDOMIZE_VARIABLES", r_val=oi_env%randomize_variables, error=error) 00661 00662 CALL section_vals_get(variable_section,explicit=explicit, n_repetition=n_var, error=error) 00663 IF (explicit) THEN 00664 ALLOCATE(oi_env%variables(1:n_var)) 00665 DO ivar=1,n_var 00666 CALL section_vals_val_get(variable_section,"VALUE",i_rep_section=ivar, & 00667 r_val=oi_env%variables(ivar)%value, error=error) 00668 CALL section_vals_val_get(variable_section,"FIXED",i_rep_section=ivar, & 00669 l_val=oi_env%variables(ivar)%fixed, error=error) 00670 CALL section_vals_val_get(variable_section,"LABEL",i_rep_section=ivar, & 00671 c_val=oi_env%variables(ivar)%label, error=error) 00672 ENDDO 00673 ENDIF 00674 00675 CALL section_vals_val_get(oi_section,"METHOD",i_val=oi_env%method,error=error) 00676 SELECT CASE(oi_env%method) 00677 CASE(opt_force_matching) 00678 fm_section => section_vals_get_subs_vals(oi_section,"FORCE_MATCHING",error=error) 00679 CALL section_vals_val_get(fm_section,"REF_TRAJ_FILE_NAME",c_val=oi_env%fm_env%ref_traj_file_name,error=error) 00680 CALL section_vals_val_get(fm_section,"REF_FORCE_FILE_NAME",c_val=oi_env%fm_env%ref_force_file_name,error=error) 00681 CALL section_vals_val_get(fm_section,"REF_CELL_FILE_NAME",c_val=oi_env%fm_env%ref_cell_file_name,error=error) 00682 CALL section_vals_val_get(fm_section,"OPTIMIZE_FILE_NAME",c_val=oi_env%fm_env%optimize_file_name,error=error) 00683 CALL section_vals_val_get(fm_section,"FRAME_START",i_val=oi_env%fm_env%frame_start,error=error) 00684 CALL section_vals_val_get(fm_section,"FRAME_STOP",i_val=oi_env%fm_env%frame_stop,error=error) 00685 CALL section_vals_val_get(fm_section,"FRAME_STRIDE",i_val=oi_env%fm_env%frame_stride,error=error) 00686 CALL section_vals_val_get(fm_section,"FRAME_COUNT",i_val=oi_env%fm_env%frame_count,error=error) 00687 00688 CALL section_vals_val_get(fm_section,"GROUP_SIZE",i_val=oi_env%fm_env%group_size,error=error) 00689 00690 CALL section_vals_val_get(fm_section,"ENERGY_WEIGHT",r_val=oi_env%fm_env%energy_weight,error=error) 00691 CALL section_vals_val_get(fm_section,"SHIFT_MM",r_val=oi_env%fm_env%shift_mm,error=error) 00692 CALL section_vals_val_get(fm_section,"SHIFT_QM",r_val=oi_env%fm_env%shift_qm,error=error) 00693 CALL section_vals_val_get(fm_section,"SHIFT_AVERAGE",l_val=oi_env%fm_env%shift_average,error=error) 00694 CASE DEFAULT 00695 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00696 END SELECT 00697 00698 CALL timestop(handle) 00699 00700 END SUBROUTINE parse_input 00701 00702 END MODULE optimize_input
1.7.3