|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00014 MODULE mc_moves 00015 USE atomic_kind_types, ONLY: get_atomic_kind 00016 USE cell_types, ONLY: cell_clone,& 00017 cell_create,& 00018 cell_release,& 00019 cell_type,& 00020 get_cell 00021 USE cp_subsys_types, ONLY: cp_subsys_get,& 00022 cp_subsys_set,& 00023 cp_subsys_type 00024 USE f77_blas 00025 USE force_env_methods, ONLY: force_env_calc_energy_force 00026 USE force_env_types, ONLY: force_env_get,& 00027 force_env_set_cell,& 00028 force_env_type,& 00029 use_fist_force 00030 USE global_types, ONLY: global_environment_type 00031 USE kinds, ONLY: default_string_length,& 00032 dp,& 00033 dp_size 00034 USE mathconstants, ONLY: pi 00035 USE mc_coordinates, ONLY: check_for_overlap,& 00036 create_discrete_array,& 00037 generate_cbmc_swap_config,& 00038 get_center_of_mass 00039 USE mc_types, ONLY: get_mc_molecule_info,& 00040 get_mc_par,& 00041 mc_ekin_type,& 00042 mc_molecule_info_type,& 00043 mc_moves_type,& 00044 mc_simpar_type 00045 USE md_run, ONLY: qs_mol_dyn 00046 USE message_passing, ONLY: mp_bcast 00047 USE mol_kind_new_list_types, ONLY: mol_kind_new_list_type 00048 USE molecule_kind_types, ONLY: bend_type,& 00049 bond_type,& 00050 get_molecule_kind,& 00051 molecule_kind_type,& 00052 torsion_type 00053 USE parallel_rng_types, ONLY: next_random_number,& 00054 rng_stream_type 00055 USE particle_list_types, ONLY: particle_list_type 00056 USE physcon, ONLY: angstrom 00057 USE termination, ONLY: stop_memory,& 00058 stop_program 00059 USE timings, ONLY: timeset,& 00060 timestop 00061 #include "cp_common_uses.h" 00062 00063 IMPLICIT NONE 00064 00065 PRIVATE 00066 00067 PRIVATE :: change_bond_angle,change_bond_length,depth_first_search,& 00068 change_dihedral 00069 00070 ! *** Global parameters *** 00071 00072 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_moves' 00073 00074 PUBLIC :: mc_conformation_change,mc_molecule_translation,& 00075 mc_molecule_rotation,mc_volume_move,mc_avbmc_move,& 00076 mc_hmc_move 00077 00078 CONTAINS 00079 00080 ! ***************************************************************************** 00094 RECURSIVE SUBROUTINE depth_first_search ( current_atom,avoid_atom,& 00095 connectivity,atom) 00096 00097 INTEGER, INTENT(IN) :: current_atom, avoid_atom 00098 INTEGER, DIMENSION(:, :), INTENT(IN) :: connectivity 00099 INTEGER, DIMENSION(:), INTENT(INOUT) :: atom 00100 00101 INTEGER :: iatom 00102 00103 DO iatom=1,6 00104 IF(connectivity(iatom,current_atom) .NE. 0) THEN 00105 IF(connectivity(iatom,current_atom) .NE. avoid_atom) THEN 00106 atom(connectivity(iatom,current_atom))=1 00107 CALL depth_first_search ( connectivity(iatom,current_atom),& 00108 current_atom,connectivity,atom) 00109 ENDIF 00110 ELSE 00111 RETURN 00112 ENDIF 00113 ENDDO 00114 00115 END SUBROUTINE depth_first_search 00116 00117 ! ***************************************************************************** 00141 SUBROUTINE mc_conformation_change ( mc_par,force_env,bias_env, moves,& 00142 move_updates,start_atom,molecule_type,box_number,& 00143 bias_energy,move_type,lreject,& 00144 rng_stream,error) 00145 00146 TYPE(mc_simpar_type), POINTER :: mc_par 00147 TYPE(force_env_type), POINTER :: force_env, bias_env 00148 TYPE(mc_moves_type), POINTER :: moves, move_updates 00149 INTEGER, INTENT(IN) :: start_atom, molecule_type, 00150 box_number 00151 REAL(KIND=dp), INTENT(INOUT) :: bias_energy 00152 CHARACTER(LEN=*), INTENT(IN) :: move_type 00153 LOGICAL, INTENT(OUT) :: lreject 00154 TYPE(rng_stream_type), POINTER :: rng_stream 00155 TYPE(cp_error_type), INTENT(INOUT) :: error 00156 00157 CHARACTER(len=*), PARAMETER :: routineN = 'mc_conformation_change', 00158 routineP = moduleN//':'//routineN 00159 00160 CHARACTER(default_string_length) :: name 00161 CHARACTER(default_string_length), 00162 DIMENSION(:), POINTER :: names 00163 INTEGER :: atom_number, end_atom, end_mol, group, handle, imol_type, 00164 imolecule, ipart, istat, jbox, molecule_number, nunits_mol, source, 00165 start_mol 00166 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits 00167 INTEGER, DIMENSION(:, :), POINTER :: nchains 00168 LOGICAL :: ionode, lbias, loverlap 00169 REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, dis_length, 00170 exp_max_val, exp_min_val, rand, value, w 00171 REAL(KIND=dp), ALLOCATABLE, 00172 DIMENSION(:, :) :: r_new, r_old 00173 TYPE(cp_subsys_type), POINTER :: subsys 00174 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00175 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds_new 00176 TYPE(molecule_kind_type), POINTER :: molecule_kind, 00177 molecule_kind_test 00178 TYPE(particle_list_type), POINTER :: particles 00179 00180 ! begin the timing of the subroutine 00181 00182 CALL timeset(routineN,handle) 00183 00184 ! nullify some pointers 00185 NULLIFY(particles,subsys,molecule_kinds_new,molecule_kind,& 00186 molecule_kind_test) 00187 00188 ! get a bunch of stuff from mc_par 00189 CALL get_mc_par(mc_par,lbias=lbias,mc_molecule_info=mc_molecule_info,& 00190 BETA=BETA,exp_max_val=exp_max_val,& 00191 exp_min_val=exp_min_val,group=group,source=source,ionode=ionode) 00192 CALL get_mc_molecule_info(mc_molecule_info,nchains=nchains,nunits=nunits,& 00193 mol_type=mol_type,names=names) 00194 00195 ! do some allocation 00196 nunits_mol=nunits(molecule_type) 00197 ALLOCATE (r_old(1:3,1:nunits_mol),STAT=istat) 00198 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00199 "r",3*nunits_mol*dp_size) 00200 ALLOCATE (r_new(1:3,1:nunits_mol),STAT=istat) 00201 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00202 "r_new",3*nunits_mol*dp_size) 00203 00204 ! find out some bounds for mol_type 00205 start_mol=1 00206 DO jbox=1,box_number-1 00207 start_mol=start_mol+SUM(nchains(:,jbox)) 00208 ENDDO 00209 end_mol=start_mol+SUM(nchains(:,box_number))-1 00210 00211 ! figure out which molecule number we are 00212 end_atom=start_atom+nunits_mol-1 00213 molecule_number=0 00214 atom_number=1 00215 DO imolecule=1,SUM(nchains(:,box_number)) 00216 IF(atom_number == start_atom) THEN 00217 molecule_number=imolecule 00218 EXIT 00219 ENDIF 00220 atom_number=atom_number+nunits(mol_type(imolecule+start_mol-1)) 00221 ENDDO 00222 IF(molecule_number == 0) CALL stop_program(routineN,moduleN,__LINE__,& 00223 'Cannot find the molecule number') 00224 00225 ! are we biasing this move? 00226 IF(lbias) THEN 00227 00228 ! grab the coordinates 00229 CALL force_env_get(bias_env,subsys=subsys,error=error) 00230 ! save the energy 00231 bias_energy_old=bias_energy 00232 00233 ELSE 00234 00235 ! grab the coordinates 00236 CALL force_env_get(force_env,subsys=subsys,error=error) 00237 ENDIF 00238 00239 ! now find the molecule type associated with this guy 00240 CALL cp_subsys_get(subsys, & 00241 particles=particles, molecule_kinds_new=molecule_kinds_new,& 00242 error=error) 00243 DO imol_type=1,SIZE(molecule_kinds_new%els(:)) 00244 molecule_kind_test => molecule_kinds_new%els(imol_type) 00245 CALL get_molecule_kind(molecule_kind_test,name=name) 00246 IF(TRIM(ADJUSTL(name)) == TRIM(ADJUSTL(names(molecule_type)))) THEN 00247 molecule_kind => molecule_kinds_new%els(imol_type) 00248 EXIT 00249 ENDIF 00250 ENDDO 00251 00252 ! save the coordinates 00253 DO ipart=start_atom,end_atom 00254 r_old(1:3,ipart-start_atom+1)=particles%els(ipart)%r(1:3) 00255 ENDDO 00256 00257 IF(.NOT. ASSOCIATED(molecule_kind)) CALL & 00258 stop_program(routineN,moduleN,__LINE__,& 00259 'Cannot find the molecule type') 00260 ! do the move 00261 IF (move_type == 'bond') THEN 00262 00263 ! record the attempt 00264 moves%bond%attempts=moves%bond%attempts+1 00265 move_updates%bond%attempts=move_updates%bond%attempts+1 00266 moves%bias_bond%attempts=moves%bias_bond%attempts+1 00267 move_updates%bias_bond%attempts=move_updates%bias_bond%attempts+1 00268 IF ( .NOT. lbias ) THEN 00269 moves%bond%qsuccesses=moves%bond%qsuccesses+1 00270 move_updates%bond%qsuccesses=& 00271 move_updates%bond%qsuccesses+1 00272 moves%bias_bond%qsuccesses=moves%bias_bond%qsuccesses+1 00273 move_updates%bias_bond%qsuccesses=& 00274 move_updates%bias_bond%qsuccesses+1 00275 ENDIF 00276 00277 ! do the move 00278 CALL change_bond_length(r_old,r_new,mc_par,molecule_type,& 00279 molecule_kind,dis_length,particles,rng_stream,error=error) 00280 00281 ELSEIF( move_type == 'angle') THEN 00282 00283 ! record the attempt 00284 moves%angle%attempts=moves%angle%attempts+1 00285 move_updates%angle%attempts=move_updates%angle%attempts+1 00286 moves%bias_angle%attempts=moves%bias_angle%attempts+1 00287 move_updates%bias_angle%attempts=move_updates%bias_angle%attempts+1 00288 IF ( .NOT. lbias ) THEN 00289 moves%angle%qsuccesses=moves%angle%qsuccesses+1 00290 move_updates%angle%qsuccesses=& 00291 move_updates%angle%qsuccesses+1 00292 moves%bias_angle%qsuccesses=moves%bias_angle%qsuccesses+1 00293 move_updates%bias_angle%qsuccesses=& 00294 move_updates%bias_angle%qsuccesses+1 00295 ENDIF 00296 00297 ! do the move 00298 CALL change_bond_angle(r_old,r_new,mc_par,molecule_type,& 00299 molecule_kind,particles,rng_stream,error=error) 00300 dis_length=1.0E0_dp 00301 ELSE 00302 ! record the attempt 00303 moves%dihedral%attempts=moves%dihedral%attempts+1 00304 move_updates%dihedral%attempts=move_updates%dihedral%attempts+1 00305 moves%bias_dihedral%attempts=moves%bias_dihedral%attempts+1 00306 move_updates%bias_dihedral%attempts=move_updates%bias_dihedral%attempts+1 00307 IF ( .NOT. lbias ) THEN 00308 moves%dihedral%qsuccesses=moves%dihedral%qsuccesses+1 00309 move_updates%dihedral%qsuccesses=& 00310 move_updates%dihedral%qsuccesses+1 00311 moves%bias_dihedral%qsuccesses=moves%bias_dihedral%qsuccesses+1 00312 move_updates%bias_dihedral%qsuccesses=& 00313 move_updates%bias_dihedral%qsuccesses+1 00314 ENDIF 00315 00316 ! do the move 00317 CALL change_dihedral(r_old,r_new,mc_par,molecule_type,& 00318 molecule_kind,particles,rng_stream,error=error) 00319 dis_length=1.0E0_dp 00320 00321 ENDIF 00322 00323 ! set the coordinates 00324 DO ipart=start_atom,end_atom 00325 particles%els(ipart)%r(1:3)=r_new(1:3,ipart-start_atom+1) 00326 ENDDO 00327 00328 ! check for overlap 00329 lreject=.FALSE. 00330 IF(lbias) THEN 00331 CALL check_for_overlap(bias_env,nchains(:,box_number),& 00332 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00333 molecule_number=molecule_number) 00334 ELSE 00335 CALL check_for_overlap(force_env,nchains(:,box_number),& 00336 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00337 molecule_number=molecule_number) 00338 IF(loverlap) lreject=.TRUE. 00339 ENDIF 00340 00341 ! if we're biasing classical, check for acceptance 00342 IF(lbias) THEN 00343 00344 ! here's where we bias the moves 00345 00346 IF(loverlap) THEN 00347 w=0.0E0_dp 00348 ELSE 00349 CALL force_env_calc_energy_force(bias_env,calc_force=.FALSE.,error=error) 00350 CALL force_env_get(bias_env,& 00351 potential_energy=bias_energy_new,error=error) 00352 ! accept or reject the move based on the Metropolis rule with a 00353 ! correction factor for the change in phase space...dis_length is 00354 ! made unitless in change_bond_length 00355 value=-BETA*(bias_energy_new-bias_energy_old) 00356 IF (value .GT. exp_max_val) THEN 00357 w=10.0_dp 00358 ELSEIF(value .LT. exp_min_val) THEN 00359 w=0.0_dp 00360 ELSE 00361 w=EXP(value)*dis_length**2 00362 ENDIF 00363 00364 ENDIF 00365 00366 IF ( w .GE. 1.0E0_dp ) THEN 00367 w=1.0E0_dp 00368 rand=0.0E0_dp 00369 ELSE 00370 IF(ionode) THEN 00371 rand=next_random_number(rng_stream,error=error) 00372 ENDIF 00373 CALL mp_bcast(rand,source,group) 00374 ENDIF 00375 00376 IF (rand .LT. w) THEN 00377 00378 ! accept the move 00379 IF (move_type == 'bond') THEN 00380 moves%bond%qsuccesses=moves%bond%qsuccesses+1 00381 move_updates%bond%successes=& 00382 move_updates%bond%successes+1 00383 moves%bias_bond%successes=moves%bias_bond%successes+1 00384 move_updates%bias_bond%successes=& 00385 move_updates%bias_bond%successes+1 00386 ELSEIF(move_type == 'angle') THEN 00387 moves%angle%qsuccesses=moves%angle%qsuccesses+1 00388 move_updates%angle%successes=& 00389 move_updates%angle%successes+1 00390 moves%bias_angle%successes=moves%bias_angle%successes+1 00391 move_updates%bias_angle%successes=& 00392 move_updates%bias_angle%successes+1 00393 ELSE 00394 moves%dihedral%qsuccesses=moves%dihedral%qsuccesses+1 00395 move_updates%dihedral%successes=& 00396 move_updates%dihedral%successes+1 00397 moves%bias_dihedral%successes=moves%bias_dihedral%successes+1 00398 move_updates%bias_dihedral%successes=& 00399 move_updates%bias_dihedral%successes+1 00400 ENDIF 00401 00402 bias_energy=bias_energy+bias_energy_new-& 00403 bias_energy_old 00404 00405 ELSE 00406 00407 ! reject the move 00408 ! restore the coordinates 00409 CALL force_env_get(bias_env,subsys=subsys,error=error) 00410 CALL cp_subsys_get(subsys,particles=particles, error=error) 00411 DO ipart=start_atom,end_atom 00412 particles%els(ipart)%r(1:3)=r_old(1:3,ipart-start_atom+1) 00413 ENDDO 00414 CALL cp_subsys_set(subsys,particles=particles,error=error) 00415 00416 ENDIF 00417 00418 ENDIF 00419 00420 ! deallocate some stuff 00421 DEALLOCATE(r_old,STAT=istat) 00422 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00423 "r") 00424 DEALLOCATE(r_new,STAT=istat) 00425 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00426 "r_new") 00427 00428 ! end the timing 00429 CALL timestop(handle) 00430 00431 END SUBROUTINE mc_conformation_change 00432 00433 ! ***************************************************************************** 00456 SUBROUTINE mc_molecule_translation( mc_par,force_env, bias_env,moves,& 00457 move_updates,start_atom,box_number,& 00458 bias_energy,molecule_type,& 00459 lreject,rng_stream,error) 00460 00461 TYPE(mc_simpar_type), POINTER :: mc_par 00462 TYPE(force_env_type), POINTER :: force_env, bias_env 00463 TYPE(mc_moves_type), POINTER :: moves, move_updates 00464 INTEGER, INTENT(IN) :: start_atom, box_number 00465 REAL(KIND=dp), INTENT(INOUT) :: bias_energy 00466 INTEGER, INTENT(IN) :: molecule_type 00467 LOGICAL, INTENT(OUT) :: lreject 00468 TYPE(rng_stream_type), POINTER :: rng_stream 00469 TYPE(cp_error_type), INTENT(INOUT) :: error 00470 00471 CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_translation', 00472 routineP = moduleN//':'//routineN 00473 00474 INTEGER :: atom_number, end_atom, end_mol, group, handle, imolecule, 00475 ipart, iparticle, istat, jbox, molecule_number, move_direction, 00476 nunits_mol, source, start_mol 00477 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot 00478 INTEGER, DIMENSION(:, :), POINTER :: nchains 00479 LOGICAL :: ionode, lbias, loverlap 00480 REAL(dp), DIMENSION(:), POINTER :: rmtrans 00481 REAL(KIND=dp) :: BETA, bias_energy_new, 00482 bias_energy_old, dis_mol, 00483 exp_max_val, exp_min_val, 00484 rand, value, w 00485 REAL(KIND=dp), ALLOCATABLE, 00486 DIMENSION(:, :) :: r_old 00487 TYPE(cp_subsys_type), POINTER :: subsys 00488 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00489 TYPE(particle_list_type), POINTER :: particles 00490 00491 ! *** Local Counters *** 00492 ! begin the timing of the subroutine 00493 00494 CALL timeset(routineN,handle) 00495 00496 ! nullify some pointers 00497 NULLIFY(particles,subsys) 00498 00499 ! get a bunch of stuff from mc_par 00500 CALL get_mc_par(mc_par,lbias=lbias,& 00501 BETA=BETA,exp_max_val=exp_max_val,& 00502 exp_min_val=exp_min_val,rmtrans=rmtrans,ionode=ionode,source=source,& 00503 group=group,mc_molecule_info=mc_molecule_info) 00504 CALL get_mc_molecule_info(mc_molecule_info,nunits_tot=nunits_tot,& 00505 nchains=nchains,nunits=nunits,mol_type=mol_type) 00506 00507 ! find out some bounds for mol_type 00508 start_mol=1 00509 DO jbox=1,box_number-1 00510 start_mol=start_mol+SUM(nchains(:,jbox)) 00511 ENDDO 00512 end_mol=start_mol+SUM(nchains(:,box_number))-1 00513 00514 ! do some allocation 00515 ALLOCATE (r_old(1:3,1:nunits_tot(box_number)),STAT=istat) 00516 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00517 "r",3*nunits_tot(box_number)*dp_size) 00518 00519 ! find the index of the last atom of this molecule, and the molecule number 00520 nunits_mol=nunits(molecule_type) 00521 end_atom=start_atom+nunits_mol-1 00522 molecule_number=0 00523 atom_number=1 00524 DO imolecule=1,SUM(nchains(:,box_number)) 00525 IF(atom_number == start_atom) THEN 00526 molecule_number=imolecule 00527 EXIT 00528 ENDIF 00529 atom_number=atom_number+nunits(mol_type(imolecule+start_mol-1)) 00530 ENDDO 00531 IF(molecule_number == 0) CALL stop_program(routineN,moduleN,__LINE__,& 00532 'Cannot find the molecule number') 00533 00534 ! are we biasing this move? 00535 IF(lbias) THEN 00536 00537 ! grab the coordinates 00538 CALL force_env_get(bias_env,subsys=subsys,error=error) 00539 CALL cp_subsys_get(subsys,particles=particles, error=error) 00540 00541 ! save the coordinates 00542 DO ipart=1,nunits_tot(box_number) 00543 r_old(1:3,ipart)=particles%els(ipart)%r(1:3) 00544 ENDDO 00545 00546 ! save the energy 00547 bias_energy_old=bias_energy 00548 00549 ELSE 00550 00551 ! grab the coordinates 00552 CALL force_env_get(force_env,subsys=subsys,error=error) 00553 CALL cp_subsys_get(subsys,particles=particles, error=error) 00554 ENDIF 00555 00556 ! record the attempt 00557 moves%trans%attempts=moves%trans%attempts+1 00558 move_updates%trans%attempts=move_updates%trans%attempts+1 00559 moves%bias_trans%attempts=moves%bias_trans%attempts+1 00560 move_updates%bias_trans%attempts=move_updates%bias_trans%attempts+1 00561 IF ( .NOT. lbias ) THEN 00562 moves%trans%qsuccesses=moves%trans%qsuccesses+1 00563 move_updates%trans%qsuccesses=move_updates%trans%qsuccesses+1 00564 moves%bias_trans%qsuccesses=moves%bias_trans%qsuccesses+1 00565 move_updates%bias_trans%qsuccesses=move_updates%bias_trans%qsuccesses+1 00566 ENDIF 00567 00568 ! move one molecule in the system 00569 00570 ! call a random number to figure out which direction we're moving 00571 IF(ionode) rand=next_random_number(rng_stream,error=error) 00572 CALL mp_bcast(rand,source,group) 00573 move_direction=AINT(3*rand)+1 00574 00575 ! call a random number to figure out how far we're moving 00576 IF(ionode) rand=next_random_number(rng_stream,error=error) 00577 CALL mp_bcast(rand,source,group) 00578 dis_mol=rmtrans(molecule_type)*(rand-0.5E0_dp)*2.0E0_dp 00579 00580 ! do the move 00581 DO iparticle=start_atom,end_atom 00582 particles%els(iparticle)%r(move_direction)=& 00583 particles%els(iparticle)%r(move_direction)+dis_mol 00584 ENDDO 00585 CALL cp_subsys_set(subsys,particles=particles,error=error) 00586 00587 ! figure out if there is any overlap...need the number of the molecule 00588 lreject=.FALSE. 00589 IF(lbias) THEN 00590 CALL check_for_overlap(bias_env,nchains(:,box_number),& 00591 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00592 molecule_number=molecule_number) 00593 ELSE 00594 CALL check_for_overlap(force_env,nchains(:,box_number),& 00595 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00596 molecule_number=molecule_number) 00597 IF(loverlap) lreject=.TRUE. 00598 ENDIF 00599 00600 ! if we're biasing with a cheaper potential, check for acceptance 00601 IF(lbias) THEN 00602 00603 ! here's where we bias the moves 00604 IF(loverlap) THEN 00605 w=0.0E0_dp 00606 ELSE 00607 CALL force_env_calc_energy_force(bias_env,calc_force=.FALSE.,error=error) 00608 CALL force_env_get(bias_env,& 00609 potential_energy=bias_energy_new,error=error) 00610 ! accept or reject the move based on the Metropolis rule 00611 value=-BETA*(bias_energy_new-bias_energy_old) 00612 IF (value .GT. exp_max_val) THEN 00613 w=10.0_dp 00614 ELSEIF(value .LT. exp_min_val) THEN 00615 w=0.0_dp 00616 ELSE 00617 w=EXP(value) 00618 ENDIF 00619 00620 ENDIF 00621 00622 IF ( w .GE. 1.0E0_dp ) THEN 00623 w=1.0E0_dp 00624 rand=0.0E0_dp 00625 ELSE 00626 IF(ionode) rand=next_random_number(rng_stream,error=error) 00627 CALL mp_bcast(rand,source,group) 00628 ENDIF 00629 00630 IF (rand .LT. w) THEN 00631 00632 ! accept the move 00633 moves%bias_trans%successes=moves%bias_trans%successes+1 00634 move_updates%bias_trans%successes=move_updates%bias_trans%successes+1 00635 moves%trans%qsuccesses=moves%trans%qsuccesses+1 00636 move_updates%trans%successes=& 00637 move_updates%trans%successes+1 00638 moves%qtrans_dis=moves%qtrans_dis+ABS(dis_mol) 00639 bias_energy=bias_energy+bias_energy_new-& 00640 bias_energy_old 00641 00642 ELSE 00643 00644 ! reject the move 00645 ! restore the coordinates 00646 CALL force_env_get(bias_env,subsys=subsys,error=error) 00647 CALL cp_subsys_get(subsys,particles=particles, error=error) 00648 DO ipart=1,nunits_tot(box_number) 00649 particles%els(ipart)%r(1:3)=r_old(1:3,ipart) 00650 ENDDO 00651 CALL cp_subsys_set(subsys,particles=particles,error=error) 00652 00653 ENDIF 00654 00655 ENDIF 00656 00657 ! deallocate some stuff 00658 DEALLOCATE(r_old,STAT=istat) 00659 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00660 "r") 00661 00662 ! end the timing 00663 CALL timestop(handle) 00664 00665 END SUBROUTINE mc_molecule_translation 00666 00667 ! ***************************************************************************** 00691 SUBROUTINE mc_molecule_rotation ( mc_par,force_env, bias_env,moves,& 00692 move_updates,box_number,& 00693 start_atom,molecule_type,bias_energy,lreject,& 00694 rng_stream,error) 00695 00696 TYPE(mc_simpar_type), POINTER :: mc_par 00697 TYPE(force_env_type), POINTER :: force_env, bias_env 00698 TYPE(mc_moves_type), POINTER :: moves, move_updates 00699 INTEGER, INTENT(IN) :: box_number, start_atom, 00700 molecule_type 00701 REAL(KIND=dp), INTENT(INOUT) :: bias_energy 00702 LOGICAL, INTENT(OUT) :: lreject 00703 TYPE(rng_stream_type), POINTER :: rng_stream 00704 TYPE(cp_error_type), INTENT(INOUT) :: error 00705 00706 CHARACTER(len=*), PARAMETER :: routineN = 'mc_molecule_rotation', 00707 routineP = moduleN//':'//routineN 00708 00709 INTEGER :: atom_number, dir, end_atom, end_mol, group, handle, ii, 00710 imolecule, ipart, istat, iunit, jbox, molecule_number, nunits_mol, 00711 source, start_mol 00712 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot 00713 INTEGER, DIMENSION(:, :), POINTER :: nchains 00714 LOGICAL :: ionode, lbias, loverlap, lx, 00715 ly 00716 REAL(dp), DIMENSION(:), POINTER :: rmrot 00717 REAL(dp), DIMENSION(:, :), POINTER :: mass 00718 REAL(KIND=dp) :: BETA, bias_energy_new, bias_energy_old, cosdg, dgamma, 00719 exp_max_val, exp_min_val, masstot, nxcm, nycm, nzcm, rand, rx, rxnew, 00720 ry, rynew, rz, rznew, sindg, value, w 00721 REAL(KIND=dp), ALLOCATABLE, 00722 DIMENSION(:, :) :: r_old 00723 TYPE(cp_subsys_type), POINTER :: subsys 00724 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00725 TYPE(particle_list_type), POINTER :: particles 00726 00727 ! begin the timing of the subroutine 00728 00729 CALL timeset(routineN,handle) 00730 00731 NULLIFY(rmrot,subsys,particles) 00732 00733 ! get a bunch of stuff from mc_par 00734 CALL get_mc_par(mc_par,lbias=lbias,& 00735 BETA=BETA,exp_max_val=exp_max_val,& 00736 exp_min_val=exp_min_val,rmrot=rmrot,mc_molecule_info=mc_molecule_info,& 00737 ionode=ionode,group=group,source=source) 00738 CALL get_mc_molecule_info(mc_molecule_info,nunits=nunits,& 00739 nunits_tot=nunits_tot,nchains=nchains,mass=mass,& 00740 mol_type=mol_type) 00741 00742 ! figure out some bounds for mol_type 00743 start_mol=1 00744 DO jbox=1,box_number-1 00745 start_mol=start_mol+SUM(nchains(:,jbox)) 00746 ENDDO 00747 end_mol=start_mol+SUM(nchains(:,box_number))-1 00748 00749 nunits_mol=nunits(molecule_type) 00750 00751 ! nullify some pointers 00752 NULLIFY(particles,subsys) 00753 00754 ! do some allocation 00755 ALLOCATE (r_old(1:3,1:nunits_tot(box_number)),STAT=istat) 00756 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00757 "r",3*nunits_tot(box_number)*dp_size) 00758 00759 ! initialize some stuff 00760 lx =.FALSE. 00761 ly =.FALSE. 00762 00763 ! determine what the final atom in the molecule is numbered, and which 00764 ! molecule number this is 00765 end_atom=start_atom+nunits_mol-1 00766 molecule_number=0 00767 atom_number=1 00768 DO imolecule=1,SUM(nchains(:,box_number)) 00769 IF(atom_number == start_atom) THEN 00770 molecule_number=imolecule 00771 EXIT 00772 ENDIF 00773 atom_number=atom_number+nunits(mol_type(imolecule+start_mol-1)) 00774 ENDDO 00775 IF(molecule_number == 0) CALL stop_program(routineN,moduleN,__LINE__,& 00776 'Cannot find the molecule number') 00777 00778 ! are we biasing this move? 00779 IF(lbias) THEN 00780 00781 ! grab the coordinates 00782 CALL force_env_get(bias_env,subsys=subsys,error=error) 00783 CALL cp_subsys_get(subsys,particles=particles, error=error) 00784 00785 ! save the coordinates 00786 DO ipart=1,nunits_tot(box_number) 00787 r_old(1:3,ipart)=particles%els(ipart)%r(1:3) 00788 ENDDO 00789 00790 ! save the energy 00791 bias_energy_old=bias_energy 00792 00793 ELSE 00794 00795 ! grab the coordinates 00796 CALL force_env_get(force_env,subsys=subsys,error=error) 00797 CALL cp_subsys_get(subsys,particles=particles, error=error) 00798 ENDIF 00799 00800 ! grab the masses 00801 masstot=SUM(mass(1:nunits(molecule_type),molecule_type)) 00802 00803 ! record the attempt 00804 moves%bias_rot%attempts=moves%bias_rot%attempts+1 00805 move_updates%bias_rot%attempts=move_updates%bias_rot%attempts+1 00806 moves%rot%attempts=moves%rot%attempts+1 00807 move_updates%rot%attempts=move_updates%rot%attempts+1 00808 IF ( .NOT. lbias ) THEN 00809 moves%rot%qsuccesses=moves%rot%qsuccesses+1 00810 move_updates%rot%qsuccesses=move_updates%rot%qsuccesses+1 00811 moves%bias_rot%qsuccesses=moves%bias_rot%qsuccesses+1 00812 move_updates%bias_rot%qsuccesses=move_updates%bias_rot%qsuccesses+1 00813 ENDIF 00814 00815 ! rotate one molecule in the system 00816 00817 ! call a random number to figure out which direction we're moving 00818 IF(ionode) rand=next_random_number(rng_stream,error=error) 00819 ! CALL RANDOM_NUMBER(rand) 00820 CALL mp_bcast(rand,source,group) 00821 dir=AINT(3*rand)+1 00822 00823 IF (dir .EQ. 1) THEN 00824 lx = .TRUE. 00825 ELSEIF (dir .EQ. 2) THEN 00826 ly = .TRUE. 00827 ENDIF 00828 00829 ! Determine new center of mass for chain i by finding the sum 00830 ! of m*r for each unit, then dividing by the total mass of the chain 00831 nxcm = 0.0E0_dp 00832 nycm = 0.0E0_dp 00833 nzcm = 0.0E0_dp 00834 DO ii = 1, nunits_mol 00835 nxcm = nxcm + particles%els(start_atom-1+ii)%r(1)* mass(ii,molecule_type) 00836 nycm = nycm + particles%els(start_atom-1+ii)%r(2)* mass(ii,molecule_type) 00837 nzcm = nzcm + particles%els(start_atom-1+ii)%r(3)* mass(ii,molecule_type) 00838 ENDDO 00839 nxcm = nxcm / masstot 00840 nycm = nycm / masstot 00841 nzcm = nzcm / masstot 00842 00843 ! call a random number to figure out how far we're moving 00844 IF(ionode) rand=next_random_number(rng_stream,error=error) 00845 CALL mp_bcast(rand,source,group) 00846 dgamma=rmrot(molecule_type)*(rand-0.5E0_dp)*2.0E0_dp 00847 00848 ! *** set up the rotation matrix *** 00849 00850 cosdg = COS( dgamma ) 00851 sindg = SIN( dgamma ) 00852 00853 IF (lx) THEN 00854 00855 ! *** ROTATE UNITS OF I AROUND X-AXIS *** 00856 00857 DO iunit = start_atom,end_atom 00858 ry = particles%els(iunit)%r(2) - nycm 00859 rz = particles%els(iunit)%r(3) - nzcm 00860 rynew = cosdg * ry - sindg * rz 00861 rznew = cosdg * rz + sindg * ry 00862 00863 particles%els(iunit)%r(2) = rynew + nycm 00864 particles%els(iunit)%r(3) = rznew + nzcm 00865 00866 ENDDO 00867 ELSEIF (ly) THEN 00868 00869 ! *** ROTATE UNITS OF I AROUND y-AXIS *** 00870 00871 DO iunit = start_atom,end_atom 00872 rx = particles%els(iunit)%r(1) - nxcm 00873 rz = particles%els(iunit)%r(3) - nzcm 00874 rxnew = cosdg * rx + sindg * rz 00875 rznew = cosdg * rz - sindg * rx 00876 00877 particles%els(iunit)%r(1) = rxnew + nxcm 00878 particles%els(iunit)%r(3) = rznew + nzcm 00879 00880 ENDDO 00881 00882 ELSE 00883 00884 ! *** ROTATE UNITS OF I AROUND z-AXIS *** 00885 00886 DO iunit = start_atom,end_atom 00887 rx = particles%els(iunit)%r(1) - nxcm 00888 ry = particles%els(iunit)%r(2) - nycm 00889 00890 rxnew = cosdg * rx - sindg * ry 00891 rynew = cosdg * ry + sindg * rx 00892 00893 particles%els(iunit)%r(1) = rxnew + nxcm 00894 particles%els(iunit)%r(2) = rynew + nycm 00895 00896 ENDDO 00897 00898 ENDIF 00899 CALL cp_subsys_set(subsys,particles=particles,error=error) 00900 00901 ! check for overlap 00902 lreject=.FALSE. 00903 IF(lbias) THEN 00904 CALL check_for_overlap(bias_env,nchains(:,box_number),& 00905 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00906 molecule_number=molecule_number) 00907 ELSE 00908 CALL check_for_overlap(force_env,nchains(:,box_number),& 00909 nunits(:),loverlap,mol_type(start_mol:end_mol),& 00910 molecule_number=molecule_number) 00911 IF(loverlap) lreject=.TRUE. 00912 ENDIF 00913 00914 ! if we're biasing classical, check for acceptance 00915 IF(lbias) THEN 00916 00917 ! here's where we bias the moves 00918 00919 IF(loverlap) THEN 00920 w=0.0E0_dp 00921 ELSE 00922 CALL force_env_calc_energy_force(bias_env,calc_force=.FALSE.,error=error) 00923 CALL force_env_get(bias_env,& 00924 potential_energy=bias_energy_new,error=error) 00925 ! accept or reject the move based on the Metropolis rule 00926 value=-BETA*(bias_energy_new-bias_energy_old) 00927 IF (value .GT. exp_max_val) THEN 00928 w=10.0_dp 00929 ELSEIF(value .LT. exp_min_val) THEN 00930 w=0.0_dp 00931 ELSE 00932 w=EXP(value) 00933 ENDIF 00934 00935 ENDIF 00936 00937 IF ( w .GE. 1.0E0_dp ) THEN 00938 w=1.0E0_dp 00939 rand=0.0E0_dp 00940 ELSE 00941 IF(ionode) rand=next_random_number(rng_stream,error=error) 00942 CALL mp_bcast(rand,source,group) 00943 ENDIF 00944 00945 IF (rand .LT. w) THEN 00946 00947 ! accept the move 00948 moves%bias_rot%successes=moves%bias_rot%successes+1 00949 move_updates%bias_rot%successes=move_updates%bias_rot%successes+1 00950 moves%rot%qsuccesses=moves%rot%qsuccesses+1 00951 move_updates%rot%successes=move_updates%rot%successes+1 00952 bias_energy=bias_energy+bias_energy_new-& 00953 bias_energy_old 00954 00955 ELSE 00956 00957 ! reject the move 00958 ! restore the coordinates 00959 CALL force_env_get(bias_env,subsys=subsys,error=error) 00960 CALL cp_subsys_get(subsys,particles=particles, error=error) 00961 DO ipart=1,nunits_tot(box_number) 00962 particles%els(ipart)%r(1:3)=r_old(1:3,ipart) 00963 ENDDO 00964 CALL cp_subsys_set(subsys,particles=particles,error=error) 00965 00966 ENDIF 00967 00968 ENDIF 00969 00970 ! deallocate some stuff 00971 DEALLOCATE(r_old,STAT=istat) 00972 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00973 "r_old") 00974 00975 ! end the timing 00976 CALL timestop(handle) 00977 00978 END SUBROUTINE mc_molecule_rotation 00979 00980 ! ***************************************************************************** 01006 SUBROUTINE mc_volume_move ( mc_par,force_env, moves,move_updates,& 01007 old_energy,box_number,& 01008 energy_check,r_old,iw,discrete_array,rng_stream,error) 01009 01010 TYPE(mc_simpar_type), POINTER :: mc_par 01011 TYPE(force_env_type), POINTER :: force_env 01012 TYPE(mc_moves_type), POINTER :: moves, move_updates 01013 REAL(KIND=dp), INTENT(INOUT) :: old_energy 01014 INTEGER, INTENT(IN) :: box_number 01015 REAL(KIND=dp), INTENT(INOUT) :: energy_check 01016 REAL(KIND=dp), DIMENSION(:, :), 01017 INTENT(INOUT) :: r_old 01018 INTEGER, INTENT(IN) :: iw 01019 INTEGER, DIMENSION(1:3, 1:2), 01020 INTENT(INOUT) :: discrete_array 01021 TYPE(rng_stream_type), POINTER :: rng_stream 01022 TYPE(cp_error_type), INTENT(INOUT) :: error 01023 01024 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_volume_move', 01025 routineP = moduleN//':'//routineN 01026 01027 CHARACTER(LEN=200) :: fft_lib 01028 CHARACTER(LEN=40) :: dat_file 01029 INTEGER :: cl, end_atom, end_mol, group, handle, iatom, idim, imolecule, 01030 iside, iside_change, istat, iunit, jbox, nunits_mol, output_unit, 01031 print_level, source, start_atom, start_mol 01032 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot 01033 INTEGER, DIMENSION(:, :), POINTER :: nchains 01034 LOGICAL :: ionode, ldiscrete, lincrease, 01035 loverlap, ltoo_small 01036 REAL(dp), DIMENSION(:, :), POINTER :: mass 01037 REAL(KIND=dp) :: BETA, discrete_step, energy_term, exp_max_val, 01038 exp_min_val, new_energy, pressure, pressure_term, rand, rcut, rmvolume, 01039 temp_var, value, vol_dis, volume_term, w 01040 REAL(KIND=dp), ALLOCATABLE, 01041 DIMENSION(:, :) :: r 01042 REAL(KIND=dp), DIMENSION(1:3) :: abc, center_of_mass, 01043 center_of_mass_new, diff, 01044 new_cell_length, 01045 old_cell_length 01046 REAL(KIND=dp), DIMENSION(1:3, 1:3) :: hmat_test 01047 TYPE(cell_type), POINTER :: cell, cell_old, cell_test 01048 TYPE(cp_logger_type), POINTER :: logger 01049 TYPE(cp_subsys_type), POINTER :: oldsys 01050 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01051 TYPE(particle_list_type), POINTER :: particles_old 01052 01053 ! begin the timing of the subroutine 01054 01055 CALL timeset(routineN,handle) 01056 01057 ! get a bunch of stuff from mc_par 01058 CALL get_mc_par(mc_par,ionode=ionode,& 01059 BETA=BETA,exp_max_val=exp_max_val,& 01060 exp_min_val=exp_min_val,source=source,group=group,& 01061 dat_file=dat_file,rmvolume=rmvolume,pressure=pressure,cl=cl,& 01062 fft_lib=fft_lib,discrete_step=discrete_step,& 01063 ldiscrete=ldiscrete,mc_molecule_info=mc_molecule_info) 01064 CALL get_mc_molecule_info(mc_molecule_info,nchains=nchains,& 01065 nunits=nunits,nunits_tot=nunits_tot,mol_type=mol_type,& 01066 mass=mass) 01067 ! figure out some bounds for mol_type 01068 start_mol=1 01069 DO jbox=1,box_number-1 01070 start_mol=start_mol+SUM(nchains(:,jbox)) 01071 ENDDO 01072 end_mol=start_mol+SUM(nchains(:,box_number))-1 01073 01074 print_level = 1 ! hack, printlevel is for print_keys 01075 01076 ! nullify some pointers 01077 NULLIFY(particles_old,cell_old,oldsys,cell_test,cell) 01078 01079 ! do some allocation 01080 ALLOCATE (r(1:3,1:nunits_tot(box_number)),STAT=istat) 01081 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01082 "r",3*nunits_tot(box_number)*dp_size) 01083 01084 ! record the attempt 01085 moves%volume%attempts=moves%volume%attempts+1 01086 move_updates%volume%attempts=move_updates%volume%attempts+1 01087 01088 ! now let's grab the cell length and particle positions 01089 CALL force_env_get(force_env,subsys=oldsys,cell=cell,& 01090 error=error) 01091 CALL get_cell(cell,abc=abc) 01092 CALL cell_create(cell_old,error=error) 01093 CALL cell_clone(cell,cell_old,error=error) 01094 CALL cp_subsys_get(oldsys,particles=particles_old, & 01095 error=error) 01096 01097 ! find the old cell length 01098 old_cell_length(1)=abc(1) 01099 old_cell_length(2)=abc(2) 01100 old_cell_length(3)=abc(3) 01101 01102 ! save the old coordiantes 01103 DO iatom=1,nunits_tot(box_number) 01104 r(1:3,iatom)=particles_old%els(iatom)%r(1:3) 01105 ENDDO 01106 01107 ! now do the move 01108 01109 ! call a random number to figure out how far we're moving 01110 IF (ionode) rand=next_random_number(rng_stream,error=error) 01111 CALL mp_bcast(rand,source,group) 01112 01113 ! find the test cell lenghts for the discrete volume move 01114 IF(ldiscrete) THEN 01115 IF(rand .LT. 0.5_dp) THEN 01116 lincrease=.TRUE. 01117 ELSE 01118 lincrease=.FALSE. 01119 ENDIF 01120 01121 new_cell_length(1:3)=old_cell_length(1:3) 01122 01123 ! if we're increasing the volume, we need to find a side we can increase 01124 IF(lincrease) THEN 01125 DO 01126 IF(ionode) rand=next_random_number(rng_stream,error=error) 01127 CALL mp_bcast(rand,source,group) 01128 iside_change=CEILING(3.0_dp*rand) 01129 IF(discrete_array(iside_change,1) .EQ. 1) THEN 01130 new_cell_length(iside_change)=& 01131 new_cell_length(iside_change)+discrete_step 01132 EXIT 01133 ENDIF 01134 ENDDO 01135 ELSE 01136 DO 01137 IF(ionode) rand=next_random_number(rng_stream,error=error) 01138 CALL mp_bcast(rand,source,group) 01139 iside_change=CEILING(3.0_dp*rand) 01140 IF(discrete_array(iside_change,2) .EQ. 1) THEN 01141 new_cell_length(iside_change)=& 01142 new_cell_length(iside_change)-discrete_step 01143 EXIT 01144 ENDIF 01145 ENDDO 01146 ENDIF 01147 vol_dis=(new_cell_length(1)*new_cell_length(2)*new_cell_length(3))& 01148 -old_cell_length(1)*old_cell_length(2)*old_cell_length(3) 01149 ELSE 01150 ! now for the not discrete volume move 01151 !!!!!!!!!!!!!!!! for E_V curves 01152 vol_dis=rmvolume*(rand-0.5E0_dp)*2.0E0_dp 01153 ! WRITE(output_unit,*) '************************ be sure to change back!',& 01154 ! old_cell_length(1),14.64_dp/angstrom 01155 ! vol_dis=-56.423592_dp/angstrom**3 01156 ! IF(old_cell_length(1) .LE. 14.64_dp/angstrom) THEN 01157 ! vol_dis=0.0_dp 01158 ! WRITE(output_unit,*) 'Found the correct box length!' 01159 ! ENDIF 01160 01161 temp_var=vol_dis+& 01162 old_cell_length(1)*old_cell_length(2)*& 01163 old_cell_length(3) 01164 01165 IF(temp_var .LE. 0.0E0_dp) THEN 01166 loverlap=.TRUE. ! cannot have a negative volume 01167 ELSE 01168 new_cell_length(1)=(temp_var)**(1.0E0_dp/3.0E0_dp) 01169 new_cell_length(2)=new_cell_length(1) 01170 new_cell_length(3)=new_cell_length(1) 01171 loverlap=.FALSE. 01172 ENDIF 01173 ENDIF 01174 CALL mp_bcast(loverlap,source,group) 01175 01176 IF(loverlap) THEN 01177 ! deallocate some stuff 01178 DEALLOCATE(r,STAT=istat) 01179 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01180 "r") 01181 logger=>cp_error_get_logger(error) 01182 output_unit=cp_logger_get_default_io_unit(logger) 01183 IF(output_unit>0) WRITE(output_unit,*) & 01184 "Volume move rejected because we tried to make too small of box.",vol_dis 01185 ! end the timing 01186 CALL timestop(handle) 01187 RETURN 01188 ENDIF 01189 01190 ! now we need to make the new cell 01191 hmat_test(:,:)=0.0e0_dp 01192 hmat_test(1,1)=new_cell_length(1) 01193 hmat_test(2,2)=new_cell_length(2) 01194 hmat_test(3,3)=new_cell_length(3) 01195 CALL cell_create(cell_test,hmat=hmat_test(:,:),periodic=cell%perd,error=error) 01196 CALL force_env_set_cell(force_env,cell_test,error=error) 01197 01198 ! now we need to scale the coordinates of all the molecules by the 01199 ! center of mass, using the minimum image (not all molecules are in 01200 ! the central box) 01201 01202 ! now we need to scale the coordinates of all the molecules by the 01203 ! center of mass 01204 end_atom=0 01205 DO imolecule=1,SUM(nchains(:,box_number)) 01206 nunits_mol=nunits(mol_type(imolecule+start_mol-1)) 01207 start_atom=end_atom+1 01208 end_atom=start_atom+nunits_mol-1 01209 ! now find the center of mass 01210 CALL get_center_of_mass(r(:,start_atom:end_atom),nunits_mol,& 01211 center_of_mass(:),mass(:,mol_type(imolecule+start_mol-1))) 01212 01213 ! scale the center of mass and determine the vector that points from the 01214 ! old COM to the new one 01215 DO iside=1,3 01216 center_of_mass_new(iside)=center_of_mass(iside)*& 01217 new_cell_length(iside)/old_cell_length(iside) 01218 ENDDO 01219 01220 DO idim=1,3 01221 diff(idim)=center_of_mass_new(idim)-center_of_mass(idim) 01222 ! now change the particle positions 01223 DO iunit=start_atom,end_atom 01224 particles_old%els(iunit)%r(idim)=& 01225 particles_old%els(iunit)%r(idim)+diff(idim) 01226 ENDDO 01227 ENDDO 01228 ENDDO 01229 01230 ! check for overlap 01231 CALL check_for_overlap(force_env,nchains(:,box_number),& 01232 nunits(:),loverlap,mol_type(start_mol:end_mol),& 01233 cell_length=new_cell_length) 01234 01235 ! figure out if we have overlap problems 01236 CALL mp_bcast(loverlap,source,group) 01237 IF(loverlap) THEN 01238 ! deallocate some stuff 01239 DEALLOCATE(r,STAT=istat) 01240 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01241 "r") 01242 01243 logger=>cp_error_get_logger(error) 01244 output_unit=cp_logger_get_default_io_unit(logger) 01245 IF(output_unit>0) WRITE(output_unit,*) & 01246 "Volume move rejected due to overlap.",vol_dis 01247 ! end the timing 01248 CALL timestop(handle) 01249 ! reset the cell and particle positions 01250 CALL force_env_set_cell(force_env,cell_old,error=error) 01251 DO iatom=1,nunits_tot(box_number) 01252 particles_old%els(iatom)%r(1:3)=r_old(1:3,iatom) 01253 ENDDO 01254 RETURN 01255 ENDIF 01256 01257 ! stop if we're trying to change a box to a boxlength smaller than rcut 01258 IF(ionode) THEN 01259 ltoo_small=.FALSE. 01260 IF(force_env%in_use .EQ. use_fist_force) THEN 01261 CALL get_mc_par(mc_par,rcut=rcut) 01262 IF(new_cell_length(1) .LT. 2.0_dp*rcut) ltoo_small=.TRUE. 01263 IF(new_cell_length(2) .LT. 2.0_dp*rcut) ltoo_small=.TRUE. 01264 IF(new_cell_length(3) .LT. 2.0_dp*rcut) ltoo_small=.TRUE. 01265 01266 IF(ltoo_small) THEN 01267 WRITE(iw,*) 'new_cell_lengths ',& 01268 new_cell_length(1:3)/angstrom 01269 WRITE(iw,*) 'rcut ',rcut/angstrom 01270 ENDIF 01271 ENDIF 01272 ENDIF 01273 CALL mp_bcast(ltoo_small,source,group) 01274 IF(ltoo_small) & 01275 CALL stop_program(routineN,moduleN,__LINE__,& 01276 "Attempted a volume move where box size got too small.") 01277 01278 ! now compute the energy 01279 CALL force_env_calc_energy_force(force_env,calc_force=.FALSE.,error=error) 01280 CALL force_env_get(force_env,& 01281 potential_energy=new_energy,error=error) 01282 01283 ! accept or reject the move 01284 ! to prevent overflows 01285 energy_term=new_energy-old_energy 01286 volume_term=-REAL(SUM(nchains(:,box_number)),dp)/BETA* 01287 LOG(new_cell_length(1)*new_cell_length(2)*new_cell_length(3)/ 01288 (old_cell_length(1)*old_cell_length(2)*old_cell_length(3))) 01289 pressure_term=pressure*vol_dis 01290 01291 value=-BETA*(energy_term+volume_term+pressure_term) 01292 IF (value .GT. exp_max_val) THEN 01293 w=10.0_dp 01294 ELSEIF(value .LT. exp_min_val) THEN 01295 w=0.0_dp 01296 ELSE 01297 w=EXP(value) 01298 ENDIF 01299 01300 !!!!!!!!!!!!!!!! for E_V curves 01301 ! w=1.0E0_dp 01302 ! w=0.0E0_dp 01303 01304 IF ( w .GE. 1.0E0_dp ) THEN 01305 w=1.0E0_dp 01306 rand=0.0E0_dp 01307 ELSE 01308 IF(ionode) rand=next_random_number(rng_stream,error=error) 01309 CALL mp_bcast(rand,source,group) 01310 ENDIF 01311 01312 IF (rand .LT. w) THEN 01313 01314 ! accept the move 01315 moves%volume%successes=moves%volume%successes+1 01316 move_updates%volume%successes=move_updates%volume%successes+1 01317 01318 ! update energies 01319 energy_check=energy_check+(new_energy-old_energy) 01320 old_energy=new_energy 01321 01322 DO iatom=1,nunits_tot(box_number) 01323 r_old(1:3,iatom)=particles_old%els(iatom)%r(1:3) 01324 ENDDO 01325 01326 ! update discrete_array if we're doing a discrete volume move 01327 IF(ldiscrete) THEN 01328 CALL create_discrete_array(new_cell_length(:),& 01329 discrete_array(:,:),discrete_step) 01330 ENDIF 01331 01332 ELSE 01333 01334 ! reset the cell and particle positions 01335 CALL force_env_set_cell(force_env,cell_old,error=error) 01336 DO iatom=1,nunits_tot(box_number) 01337 particles_old%els(iatom)%r(1:3)=r_old(1:3,iatom) 01338 ENDDO 01339 01340 ENDIF 01341 01342 ! deallocate some stuff 01343 DEALLOCATE(r,STAT=istat) 01344 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01345 "r") 01346 CALL cell_release(cell_test,error=error) 01347 CALL cell_release(cell_old,error=error) 01348 01349 ! end the timing 01350 CALL timestop(handle) 01351 01352 END SUBROUTINE mc_volume_move 01353 01354 ! ***************************************************************************** 01371 SUBROUTINE change_bond_length ( r_old,r_new,mc_par,molecule_type,molecule_kind,& 01372 dis_length,particles,rng_stream,error) 01373 01374 REAL(KIND=dp), DIMENSION(:, :), 01375 INTENT(IN) :: r_old 01376 REAL(KIND=dp), DIMENSION(:, :), 01377 INTENT(OUT) :: r_new 01378 TYPE(mc_simpar_type), POINTER :: mc_par 01379 INTEGER, INTENT(IN) :: molecule_type 01380 TYPE(molecule_kind_type), POINTER :: molecule_kind 01381 REAL(KIND=dp), INTENT(OUT) :: dis_length 01382 TYPE(particle_list_type), POINTER :: particles 01383 TYPE(rng_stream_type), POINTER :: rng_stream 01384 TYPE(cp_error_type), INTENT(inout) :: error 01385 01386 CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_length', 01387 routineP = moduleN//':'//routineN 01388 01389 INTEGER :: bond_number, group, handle, 01390 i, iatom, ibond, ipart, 01391 istat, natom, nbond, source 01392 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_b, counter 01393 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity 01394 INTEGER, DIMENSION(:), POINTER :: nunits 01395 LOGICAL :: ionode 01396 REAL(dp), DIMENSION(:), POINTER :: rmbond 01397 REAL(KIND=dp) :: atom_mass, mass_a, mass_b, 01398 new_length_a, new_length_b, 01399 old_length, rand 01400 REAL(KIND=dp), DIMENSION(1:3) :: bond_a, bond_b 01401 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 01402 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01403 01404 ! begin the timing of the subroutine 01405 01406 CALL timeset(routineN,handle) 01407 01408 NULLIFY(rmbond,mc_molecule_info) 01409 01410 ! get some stuff from mc_par 01411 CALL get_mc_par(mc_par,mc_molecule_info=mc_molecule_info,source=source,& 01412 group=group,rmbond=rmbond,ionode=ionode) 01413 CALL get_mc_molecule_info(mc_molecule_info,nunits=nunits) 01414 01415 ! copy the incoming coordinates so we can change them 01416 DO ipart=1,nunits(molecule_type) 01417 r_new(1:3,ipart)=r_old(1:3,ipart) 01418 ENDDO 01419 01420 ! pick which bond in the molecule at random 01421 IF(ionode) THEN 01422 rand=next_random_number(rng_stream,error=error) 01423 ENDIF 01424 CALL mp_bcast(rand,source,group) 01425 CALL get_molecule_kind(molecule_kind,natom=natom,nbond=nbond,& 01426 bond_list=bond_list) 01427 bond_number=CEILING(rand*REAL(nbond,dp)) 01428 01429 ALLOCATE(connection(1:natom,1:2),STAT=istat) 01430 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01431 "connection",2*natom*dp_size) 01432 ! assume at most six bonds per atom 01433 ALLOCATE(connectivity(1:6,1:natom),STAT=istat) 01434 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01435 "connectivity",6*natom*dp_size) 01436 ALLOCATE(counter(1:natom),STAT=istat) 01437 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01438 "counter",natom) 01439 ALLOCATE(atom_a(1:natom),STAT=istat) 01440 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01441 "atom_a",natom) 01442 ALLOCATE(atom_b(1:natom),STAT=istat) 01443 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01444 "atom_b",natom) 01445 connection(:,:)=0.0_dp 01446 connectivity(:,:)=0.0_dp 01447 counter(:)=0 01448 atom_a(:)=0 01449 atom_b(:)=0 01450 01451 ! now we need to find a list of atoms that each atom in this bond is connected 01452 ! to 01453 DO iatom=1,natom 01454 DO ibond=1,nbond 01455 IF(bond_list(ibond)%a == iatom) THEN 01456 counter(iatom)=counter(iatom)+1 01457 connectivity(counter(iatom),iatom)=bond_list(ibond)%b 01458 ELSEIF(bond_list(ibond)%b == iatom)THEN 01459 counter(iatom)=counter(iatom)+1 01460 connectivity(counter(iatom),iatom)=bond_list(ibond)%a 01461 ENDIF 01462 ENDDO 01463 ENDDO 01464 01465 ! now I need to do a depth first search to figure out which atoms are on atom a's 01466 ! side and which are on atom b's 01467 atom_a(:)=0 01468 atom_a(bond_list(bond_number)%a)=1 01469 CALL depth_first_search(bond_list(bond_number)%a,bond_list(bond_number)%b,& 01470 connectivity(:,:),atom_a(:)) 01471 atom_b(:)=0 01472 atom_b(bond_list(bond_number)%b)=1 01473 CALL depth_first_search(bond_list(bond_number)%b,bond_list(bond_number)%a,& 01474 connectivity(:,:),atom_b(:)) 01475 01476 ! now figure out the masses of the various sides, so we can weight how far we move each 01477 ! group of atoms 01478 mass_a=0.0_dp 01479 mass_b=0.0_dp 01480 DO iatom=1,natom 01481 CALL get_atomic_kind(particles%els(iatom)%atomic_kind,& 01482 mass=atom_mass) 01483 IF(atom_a(iatom) == 1) THEN 01484 mass_a=mass_a+atom_mass 01485 ELSE 01486 mass_b=mass_b+atom_mass 01487 ENDIF 01488 ENDDO 01489 01490 ! choose a displacement 01491 IF(ionode) rand=next_random_number(rng_stream,error=error) 01492 CALL mp_bcast(rand,source,group) 01493 01494 dis_length=rmbond(molecule_type)*2.0E0_dp*(rand-0.5E0_dp) 01495 01496 ! find the bond vector that atom a will be moving 01497 DO i=1,3 01498 bond_a(i)=r_new(i,bond_list(bond_number)%a)-& 01499 r_new(i,bond_list(bond_number)%b) 01500 bond_b(i)=-bond_a(i) 01501 ENDDO 01502 01503 ! notice we weight by the opposite masses...therefore lighter segments 01504 ! will move further 01505 old_length=SQRT(DOT_PRODUCT(bond_a,bond_a)) 01506 new_length_a=dis_length*mass_b/(mass_a+mass_b) 01507 new_length_b=dis_length*mass_a/(mass_a+mass_b) 01508 01509 DO i=1,3 01510 bond_a(i)=bond_a(i)/old_length*new_length_a 01511 bond_b(i)=bond_b(i)/old_length*new_length_b 01512 ENDDO 01513 01514 DO iatom=1,natom 01515 IF(atom_a(iatom) == 1) THEN 01516 r_new(1,iatom)=r_new(1,iatom)+bond_a(1) 01517 r_new(2,iatom)=r_new(2,iatom)+bond_a(2) 01518 r_new(3,iatom)=r_new(3,iatom)+bond_a(3) 01519 ELSE 01520 r_new(1,iatom)=r_new(1,iatom)+bond_b(1) 01521 r_new(2,iatom)=r_new(2,iatom)+bond_b(2) 01522 r_new(3,iatom)=r_new(3,iatom)+bond_b(3) 01523 ENDIF 01524 ENDDO 01525 01526 ! correct the value of dis_length for the acceptance rule 01527 dis_length=(old_length+dis_length)/old_length 01528 01529 DEALLOCATE(connection,STAT=istat) 01530 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01531 "connection") 01532 DEALLOCATE(connectivity,STAT=istat) 01533 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01534 "connectivity") 01535 DEALLOCATE(counter,STAT=istat) 01536 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01537 "counter") 01538 DEALLOCATE(atom_a,STAT=istat) 01539 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01540 "atom_a") 01541 DEALLOCATE(atom_b,STAT=istat) 01542 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01543 "atom_b") 01544 ! end the timing 01545 CALL timestop(handle) 01546 01547 END SUBROUTINE change_bond_length 01548 01549 ! ***************************************************************************** 01564 SUBROUTINE change_bond_angle ( r_old,r_new,mc_par,molecule_type,molecule_kind,& 01565 particles,rng_stream,error) 01566 01567 REAL(KIND=dp), DIMENSION(:, :), 01568 INTENT(IN) :: r_old 01569 REAL(KIND=dp), DIMENSION(:, :), 01570 INTENT(OUT) :: r_new 01571 TYPE(mc_simpar_type), POINTER :: mc_par 01572 INTEGER, INTENT(IN) :: molecule_type 01573 TYPE(molecule_kind_type), POINTER :: molecule_kind 01574 TYPE(particle_list_type), POINTER :: particles 01575 TYPE(rng_stream_type), POINTER :: rng_stream 01576 TYPE(cp_error_type), INTENT(inout) :: error 01577 01578 CHARACTER(len=*), PARAMETER :: routineN = 'change_bond_angle', 01579 routineP = moduleN//':'//routineN 01580 01581 INTEGER :: bend_number, group, handle, 01582 i, iatom, ibond, ipart, 01583 istat, natom, nbend, nbond, 01584 source 01585 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_c, counter 01586 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity 01587 INTEGER, DIMENSION(:), POINTER :: nunits 01588 LOGICAL :: ionode 01589 REAL(dp), DIMENSION(:), POINTER :: rmangle 01590 REAL(KIND=dp) :: atom_mass, bis_length, dis_angle, dis_angle_a, 01591 dis_angle_c, mass_a, mass_c, new_angle_a, new_angle_c, old_angle, 01592 old_length_a, old_length_c, rand, temp_length 01593 REAL(KIND=dp), DIMENSION(1:3) :: bisector, bond_a, bond_c, 01594 cross_prod, cross_prod_plane, 01595 temp 01596 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 01597 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 01598 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01599 01600 ! begin the timing of the subroutine 01601 01602 CALL timeset(routineN,handle) 01603 01604 NULLIFY(bend_list,bond_list,rmangle,mc_molecule_info) 01605 01606 ! get some stuff from mc_par 01607 CALL get_mc_par(mc_par,rmangle=rmangle,source=source,& 01608 group=group,ionode=ionode,mc_molecule_info=mc_molecule_info) 01609 CALL get_mc_molecule_info(mc_molecule_info,nunits=nunits) 01610 01611 ! copy the incoming coordinates so we can change them 01612 DO ipart=1,nunits(molecule_type) 01613 r_new(1:3,ipart)=r_old(1:3,ipart) 01614 ENDDO 01615 01616 ! pick which bond in the molecule at random 01617 IF(ionode) THEN 01618 rand=next_random_number(rng_stream,error=error) 01619 ENDIF 01620 CALL mp_bcast(rand,source,group) 01621 CALL get_molecule_kind(molecule_kind,natom=natom,nbend=nbend,& 01622 bend_list=bend_list,bond_list=bond_list,nbond=nbond) 01623 bend_number=CEILING(rand*REAL(nbend,dp)) 01624 01625 ALLOCATE(connection(1:natom,1:2),STAT=istat) 01626 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01627 "connection",2*natom*dp_size) 01628 ! assume at most six bonds per atom 01629 ALLOCATE(connectivity(1:6,1:natom),STAT=istat) 01630 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01631 "connectivity",6*natom*dp_size) 01632 ALLOCATE(counter(1:natom),STAT=istat) 01633 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01634 "counter",natom) 01635 ALLOCATE(atom_a(1:natom),STAT=istat) 01636 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01637 "atom_a",natom) 01638 ALLOCATE(atom_c(1:natom),STAT=istat) 01639 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01640 "atom_c",natom) 01641 connection(:,:)=0.0_dp 01642 connectivity(:,:)=0.0_dp 01643 counter(:)=0 01644 atom_a(:)=0 01645 atom_c(:)=0 01646 01647 ! now we need to find a list of atoms that each atom in this bond is connected 01648 ! to 01649 DO iatom=1,natom 01650 DO ibond=1,nbond 01651 IF(bond_list(ibond)%a == iatom) THEN 01652 counter(iatom)=counter(iatom)+1 01653 connectivity(counter(iatom),iatom)=bond_list(ibond)%b 01654 ELSEIF(bond_list(ibond)%b == iatom)THEN 01655 counter(iatom)=counter(iatom)+1 01656 connectivity(counter(iatom),iatom)=bond_list(ibond)%a 01657 ENDIF 01658 ENDDO 01659 ENDDO 01660 01661 ! now I need to do a depth first search to figure out which atoms are on atom a's 01662 ! side and which are on atom c's 01663 atom_a(:)=0 01664 atom_a(bend_list(bend_number)%a)=1 01665 CALL depth_first_search(bend_list(bend_number)%a,bend_list(bend_number)%b,& 01666 connectivity(:,:),atom_a(:)) 01667 atom_c(:)=0 01668 atom_c(bend_list(bend_number)%c)=1 01669 CALL depth_first_search(bend_list(bend_number)%c,bend_list(bend_number)%b,& 01670 connectivity(:,:),atom_c(:)) 01671 01672 ! now figure out the masses of the various sides, so we can weight how far we move each 01673 ! group of atoms 01674 mass_a=0.0_dp 01675 mass_c=0.0_dp 01676 DO iatom=1,natom 01677 CALL get_atomic_kind(particles%els(iatom)%atomic_kind,& 01678 mass=atom_mass) 01679 IF(atom_a(iatom) == 1) mass_a=mass_a+atom_mass 01680 IF(atom_c(iatom) == 1) mass_c=mass_c+atom_mass 01681 ENDDO 01682 01683 ! choose a displacement 01684 IF(ionode) rand=next_random_number(rng_stream,error=error) 01685 CALL mp_bcast(rand,source,group) 01686 01687 dis_angle=rmangle(molecule_type)*2.0E0_dp*(rand-0.5E0_dp) 01688 01689 ! need to find the A-B-C bisector 01690 01691 ! this going to be tough...we need to find the plane of the A-B-C bond and only shift 01692 ! that component for all atoms connected to A and C...otherwise we change other 01693 ! internal degrees of freedom 01694 01695 ! find the bond vectors 01696 DO i=1,3 01697 bond_a(i)=r_new(i,bend_list(bend_number)%a)-& 01698 r_new(i,bend_list(bend_number)%b) 01699 bond_c(i)=r_new(i,bend_list(bend_number)%c)-& 01700 r_new(i,bend_list(bend_number)%b) 01701 ENDDO 01702 old_length_a=SQRT(DOT_PRODUCT(bond_a,bond_a)) 01703 old_length_c=SQRT(DOT_PRODUCT(bond_c,bond_c)) 01704 old_angle=ACOS(DOT_PRODUCT(bond_a,bond_c)/(old_length_a*old_length_c)) 01705 01706 DO i=1,3 01707 bisector(i)=bond_a(i)/old_length_a+& ! not yet normalized 01708 bond_c(i)/old_length_c 01709 ENDDO 01710 bis_length=SQRT(DOT_PRODUCT(bisector,bisector)) 01711 bisector(1:3)=bisector(1:3)/bis_length 01712 01713 ! now we need to find the cross product of the B-A and B-C vectors and normalize 01714 ! it, so we have a vector that defines the bend plane 01715 cross_prod(1)=bond_a(2)*bond_c(3)-bond_a(3)*bond_c(2) 01716 cross_prod(2)=bond_a(3)*bond_c(1)-bond_a(1)*bond_c(3) 01717 cross_prod(3)=bond_a(1)*bond_c(2)-bond_a(2)*bond_c(1) 01718 cross_prod(1:3)=cross_prod(1:3)/SQRT(DOT_PRODUCT(cross_prod,cross_prod)) 01719 01720 ! we have two axis of a coordinate system...let's get the third 01721 cross_prod_plane(1)=cross_prod(2)*bisector(3)-cross_prod(3)*bisector(2) 01722 cross_prod_plane(2)=cross_prod(3)*bisector(1)-cross_prod(1)*bisector(3) 01723 cross_prod_plane(3)=cross_prod(1)*bisector(2)-cross_prod(2)*bisector(1) 01724 cross_prod_plane(1:3)=cross_prod_plane(1:3)/& 01725 SQRT(DOT_PRODUCT(cross_prod_plane,cross_prod_plane)) 01726 01727 ! now bisector is x, cross_prod_plane is the y vector (pointing towards c), 01728 ! and cross_prod is z 01729 ! shift the molecule so that atom b is at the origin 01730 DO iatom=1,natom 01731 r_new(1:3,iatom)=r_new(1:3,iatom)-& 01732 r_old(1:3,bend_list(bend_number)%b) 01733 ENDDO 01734 01735 ! figure out how much we move each side, since we're mass-weighting, by the 01736 ! opposite masses, so lighter moves farther..this angle is the angle between 01737 ! the bond vector BA or BC and the bisector 01738 dis_angle_a=dis_angle*mass_c/(mass_a+mass_c) 01739 dis_angle_c=dis_angle*mass_a/(mass_a+mass_c) 01740 01741 ! now loop through all the atoms, moving the ones that are connected to a or c 01742 DO iatom=1,natom 01743 ! subtract out the z component (perpendicular to the angle plane) 01744 temp(1:3)=r_new(1:3,iatom)-& 01745 DOT_PRODUCT(cross_prod(1:3),r_new(1:3,iatom))*& 01746 cross_prod(1:3) 01747 temp_length=SQRT(DOT_PRODUCT(temp,temp)) 01748 01749 ! we can now compute all three components of the new bond vector along the 01750 ! axis defined above 01751 IF(atom_a(iatom) == 1) THEN 01752 01753 ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector, 01754 ! as the angle computed by the dot product can't distinguish between that 01755 IF(DOT_PRODUCT(cross_prod_plane(1:3),r_new(1:3,iatom)) & 01756 .LT. 0.0_dp) THEN 01757 01758 ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is 01759 new_angle_a=ACOS(DOT_PRODUCT(bisector,temp(1:3))/& 01760 (temp_length))+dis_angle_a 01761 01762 r_new(1:3,iatom)=COS(new_angle_a)*temp_length*bisector(1:3)-& 01763 SIN(new_angle_a)*temp_length*cross_prod_plane(1:3)+& 01764 DOT_PRODUCT(cross_prod(1:3),r_new(1:3,iatom))*& 01765 cross_prod(1:3) 01766 ELSE 01767 01768 ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is 01769 new_angle_a=ACOS(DOT_PRODUCT(bisector,temp(1:3))/& 01770 (temp_length))-dis_angle_a 01771 01772 r_new(1:3,iatom)=COS(new_angle_a)*temp_length*bisector(1:3)+& 01773 SIN(new_angle_a)*temp_length*cross_prod_plane(1:3)+& 01774 DOT_PRODUCT(cross_prod(1:3),r_new(1:3,iatom))*& 01775 cross_prod(1:3) 01776 ENDIF 01777 01778 ELSEIF(atom_c(iatom) == 1) THEN 01779 01780 ! if the y-coordinate is less than zero, we need to switch the sign when we make the vector, 01781 ! as the angle computed by the dot product can't distinguish between that 01782 IF(DOT_PRODUCT(cross_prod_plane(1:3),r_new(1:3,iatom)) & 01783 .LT. 0.0_dp) THEN 01784 ! need to figure out the current iatom-B-bisector angle, so we know what the new angle is 01785 new_angle_c=ACOS(DOT_PRODUCT(bisector(1:3),temp(1:3))/& 01786 (temp_length))-dis_angle_c 01787 01788 r_new(1:3,iatom)=COS(new_angle_c)*temp_length*bisector(1:3)-& 01789 SIN(new_angle_c)*temp_length*cross_prod_plane(1:3)+& 01790 DOT_PRODUCT(cross_prod(1:3),r_new(1:3,iatom))*& 01791 cross_prod(1:3) 01792 ELSE 01793 new_angle_c=ACOS(DOT_PRODUCT(bisector(1:3),temp(1:3))/& 01794 (temp_length))+dis_angle_c 01795 01796 r_new(1:3,iatom)=COS(new_angle_c)*temp_length*bisector(1:3)+& 01797 SIN(new_angle_c)*temp_length*cross_prod_plane(1:3)+& 01798 DOT_PRODUCT(cross_prod(1:3),r_new(1:3,iatom))*& 01799 cross_prod(1:3) 01800 ENDIF 01801 ENDIF 01802 01803 ENDDO 01804 01805 DO iatom=1,natom 01806 r_new(1:3,iatom)=r_new(1:3,iatom)+& 01807 r_old(1:3,bend_list(bend_number)%b) 01808 ENDDO 01809 01810 ! deallocate some stuff 01811 DEALLOCATE(connection,STAT=istat) 01812 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01813 "connection") 01814 DEALLOCATE(connectivity,STAT=istat) 01815 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01816 "connectivity") 01817 DEALLOCATE(counter,STAT=istat) 01818 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01819 "counter") 01820 DEALLOCATE(atom_a,STAT=istat) 01821 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01822 "atom_a") 01823 DEALLOCATE(atom_c,STAT=istat) 01824 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01825 "atom_c") 01826 01827 ! end the timing 01828 CALL timestop(handle) 01829 01830 END SUBROUTINE change_bond_angle 01831 01832 ! ***************************************************************************** 01850 SUBROUTINE change_dihedral ( r_old,r_new,mc_par,molecule_type,molecule_kind,& 01851 particles,rng_stream,error) 01852 01853 REAL(KIND=dp), DIMENSION(:, :), 01854 INTENT(IN) :: r_old 01855 REAL(KIND=dp), DIMENSION(:, :), 01856 INTENT(OUT) :: r_new 01857 TYPE(mc_simpar_type), POINTER :: mc_par 01858 INTEGER, INTENT(IN) :: molecule_type 01859 TYPE(molecule_kind_type), POINTER :: molecule_kind 01860 TYPE(particle_list_type), POINTER :: particles 01861 TYPE(rng_stream_type), POINTER :: rng_stream 01862 TYPE(cp_error_type), INTENT(inout) :: error 01863 01864 CHARACTER(len=*), PARAMETER :: routineN = 'change_dihedral', 01865 routineP = moduleN//':'//routineN 01866 01867 INTEGER :: group, handle, i, iatom, 01868 ibond, ipart, istat, natom, 01869 nbond, ntorsion, source, 01870 torsion_number 01871 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_a, atom_d, counter 01872 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: connection, connectivity 01873 INTEGER, DIMENSION(:), POINTER :: nunits 01874 LOGICAL :: ionode 01875 REAL(dp), DIMENSION(:), POINTER :: rmdihedral 01876 REAL(KIND=dp) :: atom_mass, dis_angle, 01877 dis_angle_a, dis_angle_d, 01878 mass_a, mass_d, old_length_a, 01879 rand, u, v, w, x, y, z 01880 REAL(KIND=dp), DIMENSION(1:3) :: bond_a, temp 01881 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 01882 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01883 TYPE(torsion_type), DIMENSION(:), 01884 POINTER :: torsion_list 01885 01886 ! begin the timing of the subroutine 01887 01888 CALL timeset(routineN,handle) 01889 01890 NULLIFY(rmdihedral,torsion_list,bond_list,mc_molecule_info) 01891 01892 ! get some stuff from mc_par 01893 CALL get_mc_par(mc_par,rmdihedral=rmdihedral,& 01894 source=source,group=group,ionode=ionode,& 01895 mc_molecule_info=mc_molecule_info) 01896 CALL get_mc_molecule_info(mc_molecule_info,nunits=nunits) 01897 01898 ! copy the incoming coordinates so we can change them 01899 DO ipart=1,nunits(molecule_type) 01900 r_new(1:3,ipart)=r_old(1:3,ipart) 01901 ENDDO 01902 01903 ! pick which bond in the molecule at random 01904 IF(ionode) THEN 01905 rand=next_random_number(rng_stream,error=error) 01906 ! CALL RANDOM_NUMBER(rand) 01907 ENDIF 01908 CALL mp_bcast(rand,source,group) 01909 CALL get_molecule_kind(molecule_kind,natom=natom,& 01910 bond_list=bond_list,nbond=nbond,& 01911 ntorsion=ntorsion,torsion_list=torsion_list) 01912 torsion_number=CEILING(rand*REAL(ntorsion,dp)) 01913 01914 ALLOCATE(connection(1:natom,1:2),STAT=istat) 01915 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01916 "connection",2*natom*dp_size) 01917 ! assume at most six bonds per atom 01918 ALLOCATE(connectivity(1:6,1:natom),STAT=istat) 01919 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01920 "connectivity",6*natom*dp_size) 01921 ALLOCATE(counter(1:natom),STAT=istat) 01922 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01923 "counter",natom) 01924 ALLOCATE(atom_a(1:natom),STAT=istat) 01925 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01926 "atom_a",natom) 01927 ALLOCATE(atom_d(1:natom),STAT=istat) 01928 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01929 "atom_d",natom) 01930 connection(:,:)=0.0_dp 01931 connectivity(:,:)=0.0_dp 01932 counter(:)=0 01933 atom_a(:)=0 01934 atom_d(:)=0 01935 01936 ! now we need to find a list of atoms that each atom in this bond is connected 01937 ! to 01938 DO iatom=1,natom 01939 DO ibond=1,nbond 01940 IF(bond_list(ibond)%a == iatom) THEN 01941 counter(iatom)=counter(iatom)+1 01942 connectivity(counter(iatom),iatom)=bond_list(ibond)%b 01943 ELSEIF(bond_list(ibond)%b == iatom)THEN 01944 counter(iatom)=counter(iatom)+1 01945 connectivity(counter(iatom),iatom)=bond_list(ibond)%a 01946 ENDIF 01947 ENDDO 01948 ENDDO 01949 01950 ! now I need to do a depth first search to figure out which atoms are on atom 01951 ! a's side and which are on atom d's, but remember we're moving all atoms on a's 01952 ! side of b, including atoms not in a's branch 01953 atom_a(:)=0 01954 atom_a(torsion_list(torsion_number)%a)=1 01955 CALL depth_first_search(torsion_list(torsion_number)%b,& 01956 torsion_list(torsion_number)%c,connectivity(:,:),atom_a(:)) 01957 atom_d(:)=0 01958 atom_d(torsion_list(torsion_number)%d)=1 01959 CALL depth_first_search(torsion_list(torsion_number)%c,& 01960 torsion_list(torsion_number)%b,connectivity(:,:),atom_d(:)) 01961 01962 ! now figure out the masses of the various sides, so we can weight how far we 01963 ! move each group of atoms 01964 mass_a=0.0_dp 01965 mass_d=0.0_dp 01966 DO iatom=1,natom 01967 CALL get_atomic_kind(particles%els(iatom)%atomic_kind,& 01968 mass=atom_mass) 01969 IF(atom_a(iatom) == 1) mass_a=mass_a+atom_mass 01970 IF(atom_d(iatom) == 1) mass_d=mass_d+atom_mass 01971 ENDDO 01972 01973 ! choose a displacement 01974 IF(ionode) rand=next_random_number(rng_stream,error=error) 01975 CALL mp_bcast(rand,source,group) 01976 01977 dis_angle=rmdihedral(molecule_type)*2.0E0_dp*(rand-0.5E0_dp) 01978 01979 ! find the bond vectors, B-C, so we know what to rotate around 01980 DO i=1,3 01981 bond_a(i)=r_new(i,torsion_list(torsion_number)%c)-& 01982 r_new(i,torsion_list(torsion_number)%b) 01983 ENDDO 01984 old_length_a=SQRT(DOT_PRODUCT(bond_a,bond_a)) 01985 bond_a(1:3)=bond_a(1:3)/old_length_a 01986 01987 ! figure out how much we move each side, since we're mass-weighting, by the 01988 ! opposite masses, so lighter moves farther...we take the opposite sign of d 01989 ! so we're not rotating both angles in the same direction 01990 dis_angle_a=dis_angle*mass_d/(mass_a+mass_d) 01991 dis_angle_d=-dis_angle*mass_a/(mass_a+mass_d) 01992 01993 DO iatom=1,natom 01994 01995 IF(atom_a(iatom) == 1) THEN 01996 ! shift the coords so b is at the origin 01997 r_new(1:3,iatom)=r_new(1:3,iatom)-& 01998 r_new(1:3,torsion_list(torsion_number)%b) 01999 02000 ! multiply by the rotation matrix 02001 u=bond_a(1) 02002 v=bond_a(2) 02003 w=bond_a(3) 02004 x=r_new(1,iatom) 02005 y=r_new(2,iatom) 02006 z=r_new(3,iatom) 02007 temp(1)=(u*(u*x+v*y+w*z)+(x*(v**2+w**2)-u*(v*y+w*z))*COS(dis_angle_a)+& 02008 SQRT(u**2+v**2+w**2)*(v*z-w*y)*SIN(dis_angle_a))/(u**2+v**2+w**2) 02009 temp(2)=(v*(u*x+v*y+w*z)+(y*(u**2+w**2)-v*(u*x+w*z))*COS(dis_angle_a)+& 02010 SQRT(u**2+v**2+w**2)*(w*x-u*z)*SIN(dis_angle_a))/(u**2+v**2+w**2) 02011 temp(3)=(w*(u*x+v*y+w*z)+(z*(v**2+u**2)-w*(u*x+v*y))*COS(dis_angle_a)+& 02012 SQRT(u**2+v**2+w**2)*(u*y-v*x)*SIN(dis_angle_a))/(u**2+v**2+w**2) 02013 02014 ! shift back to the original position 02015 temp(1:3)=temp(1:3)+r_new(1:3,torsion_list(torsion_number)%b) 02016 r_new(1:3,iatom)=temp(1:3) 02017 02018 ELSEIF(atom_d(iatom) == 1) THEN 02019 02020 ! shift the coords so c is at the origin 02021 r_new(1:3,iatom)=r_new(1:3,iatom)-& 02022 r_new(1:3,torsion_list(torsion_number)%c) 02023 02024 ! multiply by the rotation matrix 02025 u=bond_a(1) 02026 v=bond_a(2) 02027 w=bond_a(3) 02028 x=r_new(1,iatom) 02029 y=r_new(2,iatom) 02030 z=r_new(3,iatom) 02031 temp(1)=(u*(u*x+v*y+w*z)+(x*(v**2+w**2)-u*(v*y+w*z))*COS(dis_angle_d)+& 02032 SQRT(u**2+v**2+w**2)*(v*z-w*y)*SIN(dis_angle_d))/(u**2+v**2+w**2) 02033 temp(2)=(v*(u*x+v*y+w*z)+(y*(u**2+w**2)-v*(u*x+w*z))*COS(dis_angle_d)+& 02034 SQRT(u**2+v**2+w**2)*(w*x-u*z)*SIN(dis_angle_d))/(u**2+v**2+w**2) 02035 temp(3)=(w*(u*x+v*y+w*z)+(z*(v**2+u**2)-w*(u*x+v*y))*COS(dis_angle_d)+& 02036 SQRT(u**2+v**2+w**2)*(u*y-v*x)*SIN(dis_angle_d))/(u**2+v**2+w**2) 02037 02038 ! shift back to the original position 02039 temp(1:3)=temp(1:3)+r_new(1:3,torsion_list(torsion_number)%c) 02040 r_new(1:3,iatom)=temp(1:3) 02041 ENDIF 02042 ENDDO 02043 02044 ! deallocate some stuff 02045 DEALLOCATE(connection,STAT=istat) 02046 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02047 "connection") 02048 DEALLOCATE(connectivity,STAT=istat) 02049 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02050 "connectivity") 02051 DEALLOCATE(counter,STAT=istat) 02052 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02053 "counter") 02054 DEALLOCATE(atom_a,STAT=istat) 02055 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02056 "atom_a") 02057 DEALLOCATE(atom_d,STAT=istat) 02058 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02059 "atom_d") 02060 02061 ! end the timing 02062 CALL timestop(handle) 02063 02064 END SUBROUTINE change_dihedral 02065 02066 ! ***************************************************************************** 02091 SUBROUTINE mc_avbmc_move ( mc_par,force_env,bias_env, moves,& 02092 energy_check,r_old,old_energy,start_atom_swap,& 02093 target_atom,& 02094 molecule_type,box_number,bias_energy_old,last_bias_energy,& 02095 move_type,rng_stream,error) 02096 02097 TYPE(mc_simpar_type), POINTER :: mc_par 02098 TYPE(force_env_type), POINTER :: force_env, bias_env 02099 TYPE(mc_moves_type), POINTER :: moves 02100 REAL(KIND=dp), INTENT(INOUT) :: energy_check 02101 REAL(KIND=dp), DIMENSION(:, :), 02102 INTENT(INOUT) :: r_old 02103 REAL(KIND=dp), INTENT(INOUT) :: old_energy 02104 INTEGER, INTENT(IN) :: start_atom_swap, target_atom, 02105 molecule_type, box_number 02106 REAL(KIND=dp), INTENT(INOUT) :: bias_energy_old, 02107 last_bias_energy 02108 CHARACTER(LEN=*), INTENT(IN) :: move_type 02109 TYPE(rng_stream_type), POINTER :: rng_stream 02110 TYPE(cp_error_type), INTENT(INOUT) :: error 02111 02112 CHARACTER(len=*), PARAMETER :: routineN = 'mc_avbmc_move', 02113 routineP = moduleN//':'//routineN 02114 02115 INTEGER :: end_mol, group, handle, 02116 ipart, istat, jbox, natom, 02117 nswapmoves, source, start_mol 02118 INTEGER, DIMENSION(:), POINTER :: avbmc_atom, mol_type, nunits, 02119 nunits_tot 02120 INTEGER, DIMENSION(:, :), POINTER :: nchains 02121 LOGICAL :: ionode, lbias, ldum, lin, 02122 loverlap 02123 REAL(dp), DIMENSION(:), POINTER :: avbmc_rmax, avbmc_rmin, pbias 02124 REAL(dp), DIMENSION(:, :), POINTER :: mass 02125 REAL(KIND=dp) :: BETA, bias_energy_new, del_quickstep_energy, distance, 02126 exp_max_val, exp_min_val, max_val, min_val, new_energy, prefactor, 02127 rand, rdum, volume_in, volume_out, w, weight_new, weight_old 02128 REAL(KIND=dp), ALLOCATABLE, 02129 DIMENSION(:, :) :: r_new 02130 REAL(KIND=dp), DIMENSION(1:3) :: abc, RIJ 02131 TYPE(cell_type), POINTER :: cell 02132 TYPE(cp_subsys_type), POINTER :: subsys, subsys_force 02133 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 02134 TYPE(mol_kind_new_list_type), POINTER :: molecule_kinds_new 02135 TYPE(molecule_kind_type), POINTER :: molecule_kind 02136 TYPE(particle_list_type), POINTER :: particles, particles_force 02137 02138 rdum=1.0_dp 02139 02140 ! begin the timing of the subroutine 02141 CALL timeset(routineN,handle) 02142 02143 ! get a bunch of stuff from mc_par 02144 CALL get_mc_par(mc_par,lbias=lbias,& 02145 BETA=BETA,max_val=max_val, min_val=min_val, exp_max_val=exp_max_val,& 02146 exp_min_val=exp_min_val,avbmc_atom=avbmc_atom,& 02147 avbmc_rmin=avbmc_rmin,avbmc_rmax=avbmc_rmax,& 02148 nswapmoves=nswapmoves,ionode=ionode,source=source,& 02149 group=group,pbias=pbias,mc_molecule_info=mc_molecule_info) 02150 CALL get_mc_molecule_info(mc_molecule_info,nchains=nchains,& 02151 mass=mass,nunits=nunits,nunits_tot=nunits_tot,mol_type=mol_type) 02152 ! figure out some bounds for mol_type 02153 start_mol=1 02154 DO jbox=1,box_number-1 02155 start_mol=start_mol+SUM(nchains(:,jbox)) 02156 ENDDO 02157 end_mol=start_mol+SUM(nchains(:,box_number))-1 02158 02159 ! nullify some pointers 02160 NULLIFY(particles,subsys,molecule_kinds_new,molecule_kind,& 02161 particles_force,subsys_force) 02162 02163 ! do some allocation 02164 ALLOCATE (r_new(1:3,1:nunits_tot(box_number)),STAT=istat) 02165 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02166 "r_new",3*nunits_tot(box_number)*dp_size) 02167 02168 ! now we need to grab and save coordinates, in case we reject 02169 ! are we biasing this move? 02170 IF(lbias) THEN 02171 02172 ! grab the coordinates 02173 CALL force_env_get(bias_env,cell=cell,subsys=subsys,error=error) 02174 CALL cp_subsys_get(subsys, & 02175 particles=particles, molecule_kinds_new=molecule_kinds_new,& 02176 error=error) 02177 molecule_kind => molecule_kinds_new%els(1) 02178 CALL get_molecule_kind(molecule_kind,natom=natom) 02179 CALL get_cell(cell,abc=abc) 02180 02181 ! save the energy 02182 ! bias_energy_old=bias_energy 02183 02184 ELSE 02185 02186 ! grab the coordinates 02187 CALL force_env_get(force_env,cell=cell,subsys=subsys,error=error) 02188 CALL cp_subsys_get(subsys, & 02189 particles=particles, molecule_kinds_new=molecule_kinds_new,& 02190 error=error) 02191 molecule_kind => molecule_kinds_new%els(1) 02192 CALL get_molecule_kind(molecule_kind,natom=natom) 02193 CALL get_cell(cell,abc=abc) 02194 02195 ENDIF 02196 02197 ! let's determine if the molecule to be moved is in the "in" region or the 02198 ! "out" region of the target 02199 RIJ(1)=particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(1)-& 02200 particles%els(target_atom)%r(1)-abc(1)*ANINT(& 02201 (particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(1)-& 02202 particles%els(target_atom)%r(1))/abc(1)) 02203 RIJ(2)=particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(2)-& 02204 particles%els(target_atom)%r(2)-abc(2)*ANINT(& 02205 (particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(2)-& 02206 particles%els(target_atom)%r(2))/abc(2)) 02207 RIJ(3)=particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(3)-& 02208 particles%els(target_atom)%r(3)-abc(3)*ANINT(& 02209 (particles%els(start_atom_swap+avbmc_atom(molecule_type)-1)%r(3)-& 02210 particles%els(target_atom)%r(3))/abc(3)) 02211 distance=SQRT(RIJ(1)**2+RIJ(2)**2+RIJ(3)**2) 02212 IF(distance .LE. avbmc_rmax(molecule_type) .AND. distance .GE. avbmc_rmin(molecule_type)) THEN 02213 lin=.TRUE. 02214 ELSE 02215 lin=.FALSE. 02216 ENDIF 02217 02218 ! increment the counter of the particular move we've done 02219 ! swapping into the "in" region of mol_target 02220 IF(lin) THEN 02221 IF(move_type == 'in' ) THEN 02222 moves%avbmc_inin%attempts=& 02223 moves%avbmc_inin%attempts+1 02224 ELSE 02225 moves%avbmc_inout%attempts=& 02226 moves%avbmc_inout%attempts+1 02227 ENDIF 02228 ELSE 02229 IF(move_type == 'in' ) THEN 02230 moves%avbmc_outin%attempts=& 02231 moves%avbmc_outin%attempts+1 02232 ELSE 02233 moves%avbmc_outout%attempts=& 02234 moves%avbmc_outout%attempts+1 02235 ENDIF 02236 ENDIF 02237 02238 IF(lbias) THEN 02239 02240 IF(move_type == 'in') THEN 02241 02242 ! do CBMC for the old config 02243 CALL generate_cbmc_swap_config( bias_env, BETA, max_val, min_val, exp_max_val,& 02244 exp_min_val,nswapmoves,& 02245 weight_old,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02246 mass(:,molecule_type),ldum, rdum,& 02247 bias_energy_old,ionode,.TRUE.,mol_type(start_mol:end_mol),nchains(:,box_number),& 02248 source,group,rng_stream,error,& 02249 avbmc_atom=avbmc_atom(molecule_type),& 02250 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type='out',& 02251 target_atom=target_atom) 02252 02253 ELSE 02254 02255 ! do CBMC for the old config 02256 CALL generate_cbmc_swap_config( bias_env, BETA, max_val, min_val, exp_max_val,& 02257 exp_min_val,nswapmoves,& 02258 weight_old,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02259 mass(:,molecule_type),ldum, rdum,& 02260 bias_energy_old,ionode,.TRUE.,mol_type(start_mol:end_mol),nchains(:,box_number),& 02261 source,group,rng_stream,error,& 02262 avbmc_atom=avbmc_atom(molecule_type),& 02263 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type='in',& 02264 target_atom=target_atom) 02265 02266 ENDIF 02267 02268 ! generate the new config 02269 CALL generate_cbmc_swap_config( bias_env, BETA, max_val, min_val, exp_max_val,& 02270 exp_min_val, nswapmoves,& 02271 weight_new,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02272 mass(:,molecule_type),loverlap, bias_energy_new,& 02273 bias_energy_old,ionode,.FALSE.,mol_type(start_mol:end_mol),nchains(:,box_number),& 02274 source,group,rng_stream,error,& 02275 avbmc_atom=avbmc_atom(molecule_type),& 02276 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type=move_type,& 02277 target_atom=target_atom) 02278 02279 ! the energy that comes out of the above routine is the difference...we want 02280 ! the real energy for the acceptance rule...we don't do this for the 02281 ! lbias=.false. case because it doesn't appear in the acceptance rule, and 02282 ! we compensate in case of acceptance 02283 bias_energy_new=bias_energy_new+bias_energy_old 02284 02285 ELSE 02286 02287 IF(move_type == 'in') THEN 02288 02289 ! find the weight of the old config 02290 CALL generate_cbmc_swap_config( force_env, BETA, max_val, min_val, exp_max_val,& 02291 exp_min_val,nswapmoves,& 02292 weight_old,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02293 mass(:,molecule_type),ldum,rdum,old_energy,& 02294 ionode,.TRUE.,mol_type(start_mol:end_mol),nchains(:,box_number),source,group,rng_stream,error,& 02295 avbmc_atom=avbmc_atom(molecule_type),& 02296 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type='out',& 02297 target_atom=target_atom) 02298 02299 ELSE 02300 02301 ! find the weight of the old config 02302 CALL generate_cbmc_swap_config( force_env, BETA, max_val, min_val, exp_max_val,& 02303 exp_min_val,nswapmoves,& 02304 weight_old,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02305 mass(:,molecule_type),ldum,rdum,old_energy,& 02306 ionode,.TRUE.,mol_type(start_mol:end_mol),nchains(:,box_number),source,group,rng_stream,error,& 02307 avbmc_atom=avbmc_atom(molecule_type),& 02308 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type='in',& 02309 target_atom=target_atom) 02310 02311 ENDIF 02312 02313 ! generate the new config...do this after, because it changes the force_env 02314 CALL generate_cbmc_swap_config( force_env, BETA, max_val, min_val, exp_max_val,& 02315 exp_min_val,nswapmoves,& 02316 weight_new,start_atom_swap,nunits_tot(box_number),nunits,nunits(molecule_type),& 02317 mass(:,molecule_type),loverlap,new_energy,old_energy,& 02318 ionode,.FALSE.,mol_type(start_mol:end_mol),nchains(:,box_number),source,group,rng_stream,error,& 02319 avbmc_atom=avbmc_atom(molecule_type),& 02320 rmin=avbmc_rmin(molecule_type),rmax=avbmc_rmax(molecule_type),move_type=move_type,& 02321 target_atom=target_atom) 02322 02323 ENDIF 02324 02325 IF(loverlap) THEN 02326 DEALLOCATE(r_new,STAT=istat) 02327 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02328 "r_new") 02329 02330 ! need to reset the old coordinates 02331 IF(lbias) THEN 02332 CALL force_env_get(bias_env,subsys=subsys,error=error) 02333 CALL cp_subsys_get(subsys,particles=particles,error=error) 02334 ELSE 02335 CALL force_env_get(force_env,subsys=subsys,error=error) 02336 CALL cp_subsys_get(subsys,particles=particles,error=error) 02337 ENDIF 02338 DO ipart=1,nunits_tot(box_number) 02339 particles%els(ipart)%r(1:3)=r_old(1:3,ipart) 02340 ENDDO 02341 02342 CALL timestop(handle) 02343 02344 RETURN 02345 ENDIF 02346 02347 ! if we're biasing, we need to compute the new energy with the full 02348 ! potential 02349 IF(lbias) THEN 02350 ! need to give the force_env the coords from the bias_env 02351 CALL force_env_get(force_env,subsys=subsys_force,error=error) 02352 CALL cp_subsys_get(subsys_force,particles=particles_force,error=error) 02353 CALL force_env_get(bias_env,subsys=subsys,error=error) 02354 CALL cp_subsys_get(subsys,particles=particles,error=error) 02355 DO ipart=1,nunits_tot(box_number) 02356 particles_force%els(ipart)%r(1:3)=particles%els(ipart)%r(1:3) 02357 ENDDO 02358 02359 CALL force_env_calc_energy_force(force_env,& 02360 calc_force=.FALSE.,error=error) 02361 CALL force_env_get(force_env,& 02362 potential_energy=new_energy,error=error) 02363 02364 ENDIF 02365 02366 volume_in=4.0_dp/3.0_dp*pi*(avbmc_rmax(molecule_type)**3-avbmc_rmin(molecule_type)**3) 02367 volume_out=abc(1)*abc(2)*abc(3)-volume_in 02368 02369 IF(lin .AND. move_type == 'in' .OR. & 02370 .NOT. lin .AND. move_type == 'out') THEN 02371 ! standard Metropolis rule 02372 prefactor=1.0_dp 02373 ELSEIF(.NOT. lin .AND. move_type == 'in') THEN 02374 prefactor=(1.0_dp-pbias(molecule_type))*volume_in/(pbias(molecule_type)*volume_out) 02375 ELSE 02376 prefactor=pbias(molecule_type)*volume_out/((1.0_dp-pbias(molecule_type))*volume_in) 02377 ENDIF 02378 02379 IF(lbias) THEN 02380 ! AVBMC with CBMC and a biasing potential...notice that if the biasing 02381 ! potential equals the quickstep potential, this cancels out to the 02382 ! acceptance below 02383 del_quickstep_energy=(-BETA)*(new_energy-old_energy-& 02384 (bias_energy_new-bias_energy_old)) 02385 02386 IF (del_quickstep_energy .GT. exp_max_val) THEN 02387 del_quickstep_energy=max_val 02388 ELSEIF (del_quickstep_energy .LT. exp_min_val) THEN 02389 del_quickstep_energy=0.0_dp 02390 ELSE 02391 del_quickstep_energy=EXP(del_quickstep_energy) 02392 ENDIF 02393 02394 w=prefactor*del_quickstep_energy*weight_new/weight_old 02395 02396 ELSE 02397 02398 ! AVBMC with CBMC 02399 w=prefactor*weight_new/weight_old 02400 ENDIF 02401 02402 ! check if the move is accepted 02403 IF(w .GE. 1.0E0_dp) THEN 02404 rand=0.0E0_dp 02405 ELSE 02406 IF(ionode) rand=next_random_number(rng_stream,error=error) 02407 CALL mp_bcast(rand,source,group) 02408 ENDIF 02409 02410 IF ( rand .LT. w ) THEN 02411 02412 ! accept the move 02413 02414 IF(lin) THEN 02415 IF(move_type == 'in' ) THEN 02416 moves%avbmc_inin%successes=& 02417 moves%avbmc_inin%successes+1 02418 ELSE 02419 moves%avbmc_inout%successes=& 02420 moves%avbmc_inout%successes+1 02421 ENDIF 02422 ELSE 02423 IF(move_type == 'in' ) THEN 02424 moves%avbmc_outin%successes=& 02425 moves%avbmc_outin%successes+1 02426 ELSE 02427 moves%avbmc_outout%successes=& 02428 moves%avbmc_outout%successes+1 02429 ENDIF 02430 ENDIF 02431 02432 ! we need to compensate for the fact that we take the difference in 02433 ! generate_cbmc_config to keep the exponetials small 02434 IF(.NOT. lbias) THEN 02435 new_energy=new_energy+old_energy 02436 ENDIF 02437 02438 ! update energies 02439 energy_check=energy_check+(new_energy-old_energy) 02440 old_energy=new_energy 02441 02442 ! if we're biasing the update the biasing energy 02443 IF (lbias) THEN 02444 ! need to do this outside of the routine 02445 last_bias_energy=bias_energy_new 02446 bias_energy_old=bias_energy_new 02447 ENDIF 02448 02449 ! update coordinates 02450 CALL force_env_get(force_env,subsys=subsys,error=error) 02451 CALL cp_subsys_get(subsys,particles=particles,error=error) 02452 DO ipart=1,nunits_tot(box_number) 02453 r_old(1:3,ipart)=particles%els(ipart)%r(1:3) 02454 ENDDO 02455 ELSE 02456 ! reject the move...need to restore the old coordinates 02457 IF(lbias) THEN 02458 CALL force_env_get(bias_env,subsys=subsys,error=error) 02459 CALL cp_subsys_get(subsys,particles=particles, error=error) 02460 DO ipart=1,nunits_tot(box_number) 02461 particles%els(ipart)%r(1:3)=r_old(1:3,ipart) 02462 ENDDO 02463 CALL cp_subsys_set(subsys,particles=particles,error=error) 02464 ENDIF 02465 CALL force_env_get(force_env,subsys=subsys,error=error) 02466 CALL cp_subsys_get(subsys,particles=particles, error=error) 02467 DO ipart=1,nunits_tot(box_number) 02468 particles%els(ipart)%r(1:3)=r_old(1:3,ipart) 02469 ENDDO 02470 CALL cp_subsys_set(subsys,particles=particles,error=error) 02471 02472 ENDIF 02473 02474 ! deallocate some stuff 02475 DEALLOCATE(r_new,STAT=istat) 02476 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02477 "r_new") 02478 ! end the timing 02479 CALL timestop(handle) 02480 02481 END SUBROUTINE mc_avbmc_move 02482 02483 ! ***************************************************************************** 02507 SUBROUTINE mc_hmc_move ( mc_par,force_env, globenv, moves,move_updates,& 02508 old_energy,box_number,& 02509 energy_check,r_old,iw,rng_stream,error) 02510 02511 TYPE(mc_simpar_type), POINTER :: mc_par 02512 TYPE(force_env_type), POINTER :: force_env 02513 TYPE(global_environment_type), POINTER :: globenv 02514 TYPE(mc_moves_type), POINTER :: moves, move_updates 02515 REAL(KIND=dp), INTENT(INOUT) :: old_energy 02516 INTEGER, INTENT(IN) :: box_number 02517 REAL(KIND=dp), INTENT(INOUT) :: energy_check 02518 REAL(KIND=dp), DIMENSION(:, :), 02519 INTENT(INOUT) :: r_old 02520 INTEGER, INTENT(IN) :: iw 02521 TYPE(rng_stream_type), POINTER :: rng_stream 02522 TYPE(cp_error_type), INTENT(INOUT) :: error 02523 02524 CHARACTER(LEN=*), PARAMETER :: routineN = 'mc_hmc_move', 02525 routineP = moduleN//':'//routineN 02526 02527 INTEGER :: group, handle, iatom, istat, 02528 source 02529 INTEGER, DIMENSION(:), POINTER :: nunits_tot 02530 LOGICAL :: ionode 02531 REAL(KIND=dp) :: BETA, energy_term, 02532 exp_max_val, exp_min_val, 02533 new_energy, rand, value, w 02534 REAL(KIND=dp), ALLOCATABLE, 02535 DIMENSION(:, :) :: r 02536 TYPE(cp_subsys_type), POINTER :: oldsys 02537 TYPE(mc_ekin_type), POINTER :: hmc_ekin 02538 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 02539 TYPE(particle_list_type), POINTER :: particles_old 02540 02541 ! begin the timing of the subroutine 02542 02543 CALL timeset(routineN,handle) 02544 02545 ! get a bunch of stuff from mc_par 02546 CALL get_mc_par(mc_par,ionode=ionode,& 02547 BETA=BETA,exp_max_val=exp_max_val,& 02548 exp_min_val=exp_min_val,source=source,group=group,& 02549 mc_molecule_info=mc_molecule_info) 02550 CALL get_mc_molecule_info(mc_molecule_info,nunits_tot=nunits_tot) 02551 02552 ! nullify some pointers 02553 NULLIFY(particles_old,oldsys,hmc_ekin) 02554 02555 ! do some allocation 02556 ALLOCATE (r(1:3,1:nunits_tot(box_number)),STAT=istat) 02557 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02558 "r",3*nunits_tot(box_number)*dp_size) 02559 ALLOCATE (hmc_ekin,STAT=istat) 02560 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02561 "hmc_ekin") 02562 02563 ! record the attempt 02564 moves%hmc%attempts=moves%hmc%attempts+1 02565 move_updates%hmc%attempts=move_updates%hmc%attempts+1 02566 02567 ! now let's grab the particle positions 02568 CALL force_env_get(force_env,subsys=oldsys,& 02569 error=error) 02570 CALL cp_subsys_get(oldsys,particles=particles_old, & 02571 error=error) 02572 02573 ! save the old coordinates 02574 DO iatom=1,nunits_tot(box_number) 02575 r(1:3,iatom)=particles_old%els(iatom)%r(1:3) 02576 ENDDO 02577 02578 ! now run the MD simulation 02579 CALL qs_mol_dyn(force_env,globenv,error=error,hmc_ekin=hmc_ekin) 02580 02581 02582 ! get the energy 02583 CALL force_env_get(force_env,& 02584 potential_energy=new_energy,error=error) 02585 02586 ! accept or reject the move 02587 ! to prevent overflows 02588 energy_term=new_energy+hmc_ekin%final_ekin-old_energy-hmc_ekin%initial_ekin 02589 02590 value=-BETA*(energy_term) 02591 IF (value .GT. exp_max_val) THEN 02592 w=10.0_dp 02593 ELSEIF(value .LT. exp_min_val) THEN 02594 w=0.0_dp 02595 ELSE 02596 w=EXP(value) 02597 ENDIF 02598 02599 IF ( w .GE. 1.0E0_dp ) THEN 02600 w=1.0E0_dp 02601 rand=0.0E0_dp 02602 ELSE 02603 IF(ionode) rand=next_random_number(rng_stream,error=error) 02604 CALL mp_bcast(rand,source,group) 02605 ENDIF 02606 02607 IF (rand .LT. w) THEN 02608 02609 ! accept the move 02610 moves%hmc%successes=moves%hmc%successes+1 02611 move_updates%hmc%successes=move_updates%hmc%successes+1 02612 02613 ! update energies 02614 energy_check=energy_check+(new_energy-old_energy) 02615 old_energy=new_energy 02616 02617 DO iatom=1,nunits_tot(box_number) 02618 r_old(1:3,iatom)=particles_old%els(iatom)%r(1:3) 02619 ENDDO 02620 02621 ELSE 02622 02623 ! reset the cell and particle positions 02624 DO iatom=1,nunits_tot(box_number) 02625 particles_old%els(iatom)%r(1:3)=r_old(1:3,iatom) 02626 ENDDO 02627 02628 ENDIF 02629 02630 ! deallocate some stuff 02631 DEALLOCATE(r,STAT=istat) 02632 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02633 "r") 02634 DEALLOCATE(hmc_ekin,STAT=istat) 02635 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 02636 "hmc_ekin") 02637 02638 ! end the timing 02639 CALL timestop(handle) 02640 02641 END SUBROUTINE mc_hmc_move 02642 END MODULE mc_moves 02643
1.7.3