|
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 ! ***************************************************************************** 00016 MODULE mc_ge_moves 00017 USE cell_types, ONLY: cell_clone,& 00018 cell_create,& 00019 cell_p_type,& 00020 cell_release,& 00021 cell_type,& 00022 get_cell 00023 USE cp_para_types, ONLY: cp_para_env_type 00024 USE cp_subsys_types, ONLY: cp_subsys_get,& 00025 cp_subsys_p_type,& 00026 cp_subsys_set,& 00027 cp_subsys_type 00028 USE f77_blas 00029 USE force_env_methods, ONLY: force_env_calc_energy_force 00030 USE force_env_types, ONLY: force_env_get,& 00031 force_env_p_type,& 00032 force_env_release,& 00033 force_env_set_cell 00034 USE input_constants, ONLY: dump_xmol 00035 USE input_section_types, ONLY: section_vals_release 00036 USE kinds, ONLY: default_string_length,& 00037 dp,& 00038 dp_size 00039 USE mc_control, ONLY: mc_create_force_env 00040 USE mc_coordinates, ONLY: check_for_overlap,& 00041 generate_cbmc_swap_config,& 00042 get_center_of_mass 00043 USE mc_misc, ONLY: mc_make_dat_file_new 00044 USE mc_move_control, ONLY: move_q_reinit,& 00045 q_move_accept 00046 USE mc_types, ONLY: & 00047 get_mc_molecule_info, get_mc_par, mc_determine_molecule_info, & 00048 mc_input_file_type, mc_molecule_info_destroy, mc_molecule_info_type, & 00049 mc_moves_p_type, mc_simulation_parameters_p_type, set_mc_par 00050 USE message_passing, ONLY: mp_bcast 00051 USE parallel_rng_types, ONLY: next_random_number,& 00052 rng_stream_type 00053 USE particle_list_types, ONLY: particle_list_p_type,& 00054 particle_list_type 00055 USE particle_types, ONLY: write_particle_coordinates 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 ! *** Global parameters *** 00068 00069 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_ge_moves' 00070 00071 PUBLIC :: mc_ge_volume_move,mc_ge_swap_move,& 00072 mc_quickstep_move 00073 00074 CONTAINS 00075 00076 ! ***************************************************************************** 00107 SUBROUTINE mc_Quickstep_move(mc_par,force_env,bias_env,moves,& 00108 lreject,move_updates,energy_check,r_old,& 00109 nnstep,old_energy,bias_energy_new,last_bias_energy,& 00110 nboxes,box_flag,subsys,particles,rng_stream,& 00111 unit_conv,error) 00112 00113 TYPE(mc_simulation_parameters_p_type), 00114 DIMENSION(:), POINTER :: mc_par 00115 TYPE(force_env_p_type), DIMENSION(:), 00116 POINTER :: force_env, bias_env 00117 TYPE(mc_moves_p_type), DIMENSION(:, :), 00118 POINTER :: moves 00119 LOGICAL, INTENT(IN) :: lreject 00120 TYPE(mc_moves_p_type), DIMENSION(:, :), 00121 POINTER :: move_updates 00122 REAL(KIND=dp), DIMENSION(:), 00123 INTENT(INOUT) :: energy_check 00124 REAL(KIND=dp), DIMENSION(:, :, :), 00125 INTENT(INOUT) :: r_old 00126 INTEGER, INTENT(IN) :: nnstep 00127 REAL(KIND=dp), DIMENSION(:), 00128 INTENT(INOUT) :: old_energy, bias_energy_new, 00129 last_bias_energy 00130 INTEGER, INTENT(IN) :: nboxes 00131 INTEGER, DIMENSION(:), INTENT(IN) :: box_flag 00132 TYPE(cp_subsys_p_type), DIMENSION(:), 00133 POINTER :: subsys 00134 TYPE(particle_list_p_type), 00135 DIMENSION(:), POINTER :: particles 00136 TYPE(rng_stream_type), POINTER :: rng_stream 00137 REAL(KIND=dp), INTENT(IN) :: unit_conv 00138 TYPE(cp_error_type), INTENT(INOUT) :: error 00139 00140 CHARACTER(len=*), PARAMETER :: routineN = 'mc_Quickstep_move', 00141 routineP = moduleN//':'//routineN 00142 00143 INTEGER :: end_mol, group, handle, ibox, iparticle, iprint, istat, itype, 00144 jbox, nmol_types, source, start_mol 00145 INTEGER, DIMENSION(:, :), POINTER :: nchains 00146 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot 00147 INTEGER, DIMENSION(1:nboxes) :: diff 00148 LOGICAL :: ionode, lbias, loverlap 00149 REAL(KIND=dp) :: BETA, energies, rand, w 00150 REAL(KIND=dp), DIMENSION(1:nboxes) :: bias_energy_old, new_energy 00151 TYPE(cp_subsys_p_type), DIMENSION(:), 00152 POINTER :: subsys_bias 00153 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 00154 TYPE(particle_list_p_type), 00155 DIMENSION(:), POINTER :: particles_bias 00156 00157 ! begin the timing of the subroutine 00158 00159 CALL timeset(routineN,handle) 00160 00161 NULLIFY(subsys_bias,particles_bias) 00162 00163 ! get a bunch of data from mc_par 00164 CALL get_mc_par(mc_par(1)%mc_par,ionode=ionode,lbias=lbias,& 00165 BETA=BETA,diff=diff(1),source=source,group=group,& 00166 iprint=iprint,& 00167 mc_molecule_info=mc_molecule_info) 00168 CALL get_mc_molecule_info(mc_molecule_info,nmol_types=nmol_types,& 00169 nchains=nchains,nunits_tot=nunits_tot,nunits=nunits,mol_type=mol_type) 00170 00171 IF(nboxes .GT. 1) THEN 00172 DO ibox=2,nboxes 00173 CALL get_mc_par(mc_par(ibox)%mc_par,diff=diff(ibox)) 00174 ENDDO 00175 ENDIF 00176 00177 ! allocate some stuff 00178 ALLOCATE(subsys_bias(1:nboxes),STAT=istat) 00179 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00180 "subsys_bias",nboxes) 00181 ALLOCATE(particles_bias(1:nboxes),STAT=istat) 00182 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00183 "particles_bias",nboxes) 00184 00185 ! record the attempt...we really only care about molecule type 1 and box 00186 ! type 1, since the acceptance will be identical for all boxes and molecules 00187 moves(1,1)%moves%Quickstep%attempts=& 00188 moves(1,1)%moves%Quickstep%attempts+1 00189 00190 ! grab the coordinates for the force_env 00191 DO ibox=1,nboxes 00192 CALL force_env_get(force_env(ibox)%force_env,& 00193 subsys=subsys(ibox)%subsys,error=error) 00194 CALL cp_subsys_get(subsys(ibox)%subsys, & 00195 particles=particles(ibox)%list, error=error) 00196 ENDDO 00197 00198 ! calculate the new energy of the system...if we're biasing, 00199 ! force_env hasn't changed but bias_env has 00200 DO ibox=1,nboxes 00201 IF(box_flag(ibox) == 1) THEN 00202 IF(lbias) THEN 00203 ! grab the coords from bias_env and put them into force_env 00204 CALL force_env_get(bias_env(ibox)%force_env,& 00205 subsys=subsys_bias(ibox)%subsys,error=error) 00206 CALL cp_subsys_get(subsys_bias(ibox)%subsys, & 00207 particles=particles_bias(ibox)%list, error=error) 00208 00209 DO iparticle=1,nunits_tot(ibox) 00210 particles(ibox)%list%els(iparticle)%r(1:3)=& 00211 particles_bias(ibox)%list%els(iparticle)%r(1:3) 00212 ENDDO 00213 00214 CALL force_env_calc_energy_force(force_env(ibox)%force_env,& 00215 calc_force=.FALSE.,error=error) 00216 CALL force_env_get(force_env(ibox)%force_env,& 00217 potential_energy=new_energy(ibox),error=error) 00218 ELSE 00219 IF( .NOT. lreject) THEN 00220 CALL force_env_calc_energy_force(force_env(ibox)%force_env,& 00221 calc_force=.FALSE.,error=error) 00222 CALL force_env_get(force_env(ibox)%force_env,& 00223 potential_energy=new_energy(ibox),error=error) 00224 ENDIF 00225 ENDIF 00226 ELSE 00227 new_energy(ibox)=old_energy(ibox) 00228 ENDIF 00229 00230 ENDDO 00231 00232 ! accept or reject the move based on Metropolis or the Iftimie rule 00233 IF (ionode) THEN 00234 00235 ! write them out in case something bad happens 00236 IF(MOD(nnstep,iprint) == 0) THEN 00237 DO ibox=1,nboxes 00238 IF(SUM(nchains(:,ibox)) == 0) THEN 00239 WRITE(diff(ibox),*) nnstep 00240 WRITE(diff(ibox),*) nchains(:,ibox) 00241 ELSE 00242 WRITE(diff(ibox),*) nnstep 00243 CALL write_particle_coordinates(& 00244 particles(ibox)%list%els,& 00245 diff(ibox),dump_xmol,'POS','TRIAL',& 00246 unit_conv=unit_conv,error=error) 00247 ENDIF 00248 ENDDO 00249 ENDIF 00250 ENDIF 00251 00252 IF(.NOT. lreject) THEN 00253 IF (lbias) THEN 00254 00255 DO ibox=1,nboxes 00256 ! look for overlap 00257 IF(SUM(nchains(:,ibox)) .NE. 0) THEN 00258 ! find the molecule bounds 00259 start_mol=1 00260 DO jbox=1,ibox-1 00261 start_mol=start_mol+SUM(nchains(:,jbox)) 00262 ENDDO 00263 end_mol=start_mol+SUM(nchains(:,ibox))-1 00264 CALL check_for_overlap(bias_env(ibox)%force_env,& 00265 nchains(:,ibox),nunits(:),loverlap,mol_type(start_mol:end_mol)) 00266 IF(loverlap) CALL stop_program(routineN,& 00267 moduleN,__LINE__,& 00268 'Quickstep move found an overlap in the old config') 00269 ENDIF 00270 bias_energy_old(ibox)=last_bias_energy(ibox) 00271 ENDDO 00272 00273 energies=-BETA*((SUM(new_energy(:))-SUM(bias_energy_new(:)))& 00274 -(SUM(old_energy(:))-SUM(bias_energy_old(:)))) 00275 00276 ! used to prevent over and underflows 00277 IF(energies .GE. -1.0E-8) THEN 00278 w=1.0_dp 00279 ELSEIF(energies .LE. -500.0_dp) THEN 00280 w=0.0_dp 00281 ELSE 00282 w=EXP(energies) 00283 ENDIF 00284 00285 IF(ionode) THEN 00286 DO ibox=1,nboxes 00287 WRITE(diff(ibox),*) nnstep,new_energy(ibox)-& 00288 old_energy(ibox),& 00289 bias_energy_new(ibox)-bias_energy_old(ibox) 00290 ENDDO 00291 ENDIF 00292 ELSE 00293 energies=-BETA*(SUM(new_energy(:))-SUM(old_energy(:))) 00294 ! used to prevent over and underflows 00295 IF(energies .GE. 0.0_dp) THEN 00296 w=1.0_dp 00297 ELSEIF(energies .LE. -500.0_dp) THEN 00298 w=0.0_dp 00299 ELSE 00300 w=EXP(energies) 00301 ENDIF 00302 ENDIF 00303 ELSE 00304 w=0.0E0_dp 00305 ENDIF 00306 IF ( w .GE. 1.0E0_dp ) THEN 00307 w=1.0E0_dp 00308 rand=0.0E0_dp 00309 ELSE 00310 IF(ionode) rand=next_random_number(rng_stream,error=error) 00311 CALL mp_bcast(rand,source,group) 00312 ENDIF 00313 00314 IF (rand .LT. w) THEN 00315 00316 ! accept the move 00317 moves(1,1)%moves%Quickstep%successes=& 00318 moves(1,1)%moves%Quickstep%successes+1 00319 00320 DO ibox=1,nboxes 00321 ! remember what kind of move we did for lbias=.false. 00322 IF(.NOT. lbias) THEN 00323 DO itype=1,nmol_types 00324 CALL q_move_accept(moves(itype,ibox)%moves,.TRUE.) 00325 CALL q_move_accept(move_updates(itype,ibox)%moves,.TRUE.) 00326 00327 ! reset the counters 00328 CALL move_q_reinit(moves(itype,ibox)%moves, .TRUE.) 00329 CALL move_q_reinit(move_updates(itype,ibox)%moves, .TRUE.) 00330 ENDDO 00331 ENDIF 00332 00333 DO itype=1,nmol_types 00334 ! we need to record all accepted moves since last Quickstep calculation 00335 CALL q_move_accept(moves(itype,ibox)%moves,.FALSE.) 00336 CALL q_move_accept(move_updates(itype,ibox)%moves,.FALSE.) 00337 00338 ! reset the counters 00339 CALL move_q_reinit(moves(itype,ibox)%moves, .FALSE.) 00340 CALL move_q_reinit(move_updates(itype,ibox)%moves, .FALSE.) 00341 ENDDO 00342 00343 ! update energies 00344 energy_check(ibox)=energy_check(ibox)+& 00345 (new_energy(ibox)-old_energy(ibox)) 00346 old_energy(ibox)=new_energy(ibox) 00347 00348 ENDDO 00349 00350 IF ( lbias) THEN 00351 DO ibox=1,nboxes 00352 last_bias_energy(ibox)=bias_energy_new(ibox) 00353 ENDDO 00354 ENDIF 00355 00356 ! update coordinates 00357 DO ibox=1,nboxes 00358 IF(nunits_tot(ibox) .NE. 0) THEN 00359 DO iparticle=1,nunits_tot(ibox) 00360 r_old(1:3,iparticle,ibox)=& 00361 particles(ibox)%list%els(iparticle)%r(1:3) 00362 ENDDO 00363 ENDIF 00364 ENDDO 00365 00366 ELSE 00367 00368 ! reject the move 00369 DO ibox=1,nboxes 00370 DO itype=1,nmol_types 00371 CALL move_q_reinit(moves(itype,ibox)%moves,.FALSE.) 00372 CALL move_q_reinit(move_updates(itype,ibox)%moves,.FALSE.) 00373 IF(.NOT. lbias) THEN 00374 ! reset the counters 00375 CALL move_q_reinit(moves(itype,ibox)%moves, .TRUE.) 00376 CALL move_q_reinit(move_updates(itype,ibox)%moves, .TRUE.) 00377 ENDIF 00378 ENDDO 00379 00380 ENDDO 00381 00382 IF ( .NOT. ionode) r_old(:,:,:) = 0.0E0_dp 00383 00384 ! coodinates changed, so we need to broadcast those, even for the lbias 00385 ! case since bias_env needs to have the same coords as force_env 00386 CALL mp_bcast(r_old,source,group) 00387 00388 DO ibox=1,nboxes 00389 DO iparticle=1,nunits_tot(ibox) 00390 particles(ibox)%list%els(iparticle)%r(1:3)=& 00391 r_old(1:3,iparticle,ibox) 00392 IF(lbias .AND. box_flag(ibox) == 1) & 00393 particles_bias(ibox)%list%els(iparticle)%r(1:3)=& 00394 r_old(1:3,iparticle,ibox) 00395 ENDDO 00396 ENDDO 00397 00398 ! need to reset the energies of the biasing potential 00399 IF ( lbias) THEN 00400 DO ibox=1,nboxes 00401 bias_energy_new(ibox)=last_bias_energy(ibox) 00402 ENDDO 00403 ENDIF 00404 00405 ENDIF 00406 00407 ! make sure the coordinates are transferred 00408 DO ibox=1,nboxes 00409 CALL cp_subsys_set(subsys(ibox)%subsys,& 00410 particles=particles(ibox)%list,error=error) 00411 IF(lbias .AND. box_flag(ibox) == 1) & 00412 CALL cp_subsys_set(subsys_bias(ibox)%subsys,& 00413 particles=particles_bias(ibox)%list,error=error) 00414 ENDDO 00415 00416 ! deallocate some stuff 00417 DEALLOCATE(subsys_bias,STAT=istat) 00418 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00419 __LINE__,"subsys_bias") 00420 DEALLOCATE(particles_bias,STAT=istat) 00421 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00422 __LINE__,"particles_bias") 00423 00424 ! end the timing 00425 CALL timestop(handle) 00426 00427 END SUBROUTINE mc_Quickstep_move 00428 00429 ! ***************************************************************************** 00452 SUBROUTINE mc_ge_swap_move(mc_par,force_env,bias_env,moves,& 00453 energy_check,r_old,old_energy,& 00454 para_env,bias_energy_old,last_bias_energy,& 00455 rng_stream,error) 00456 00457 TYPE(mc_simulation_parameters_p_type), 00458 DIMENSION(:), POINTER :: mc_par 00459 TYPE(force_env_p_type), DIMENSION(:), 00460 POINTER :: force_env, bias_env 00461 TYPE(mc_moves_p_type), DIMENSION(:, :), 00462 POINTER :: moves 00463 REAL(KIND=dp), DIMENSION(1:2), 00464 INTENT(INOUT) :: energy_check 00465 REAL(KIND=dp), DIMENSION(:, :, :), 00466 INTENT(INOUT) :: r_old 00467 REAL(KIND=dp), DIMENSION(1:2), 00468 INTENT(INOUT) :: old_energy 00469 TYPE(cp_para_env_type), POINTER :: para_env 00470 REAL(KIND=dp), DIMENSION(1:2), 00471 INTENT(INOUT) :: bias_energy_old, 00472 last_bias_energy 00473 TYPE(rng_stream_type), POINTER :: rng_stream 00474 TYPE(cp_error_type), INTENT(INOUT) :: error 00475 00476 CHARACTER(len=*), PARAMETER :: routineN = 'mc_ge_swap_move', 00477 routineP = moduleN//':'//routineN 00478 00479 CHARACTER(default_string_length), 00480 ALLOCATABLE, DIMENSION(:) :: atom_names_insert, 00481 atom_names_remove 00482 CHARACTER(default_string_length), 00483 DIMENSION(:, :), POINTER :: atom_names 00484 CHARACTER(LEN=200) :: fft_lib 00485 CHARACTER(LEN=40), DIMENSION(1:2) :: dat_file 00486 INTEGER :: end_mol, group, handle, iatom, ibox, idim, iiatom, imolecule, 00487 ins_atoms, insert_box, ipart, istat, itype, jbox, molecule_type, 00488 nmol_types, nswapmoves, print_level, rem_atoms, remove_box, source, 00489 start_atom_ins, start_atom_rem, start_mol 00490 INTEGER, DIMENSION(:), POINTER :: mol_type, mol_type_test, 00491 nunits, nunits_tot 00492 INTEGER, DIMENSION(:, :), POINTER :: nchains, nchains_test 00493 LOGICAL :: ionode, lbias, loverlap, 00494 loverlap_ins, loverlap_rem 00495 REAL(dp), DIMENSION(:), POINTER :: eta_insert, eta_remove, 00496 pmswap_mol 00497 REAL(dp), DIMENSION(:, :), POINTER :: insert_coords, remove_coords 00498 REAL(KIND=dp) :: BETA, del_quickstep_energy, exp_max_val, exp_min_val, 00499 max_val, min_val, prefactor, rand, rdum, w, weight_new, weight_old 00500 REAL(KIND=dp), ALLOCATABLE, 00501 DIMENSION(:, :) :: cbmc_energies, r_cbmc, 00502 r_insert_mol 00503 REAL(KIND=dp), DIMENSION(1:2) :: bias_energy_new, new_energy 00504 REAL(KIND=dp), DIMENSION(1:3) :: abc_insert, abc_remove, 00505 center_of_mass, 00506 displace_molecule, pos_insert 00507 REAL(KIND=dp), DIMENSION(:, :), POINTER :: mass 00508 TYPE(cell_type), POINTER :: cell_insert, cell_remove 00509 TYPE(cp_subsys_p_type), DIMENSION(:), 00510 POINTER :: oldsys 00511 TYPE(cp_subsys_type), POINTER :: insert_sys, remove_sys 00512 TYPE(force_env_p_type), DIMENSION(:), 00513 POINTER :: test_env, test_env_bias 00514 TYPE(mc_input_file_type), POINTER :: mc_bias_file, mc_input_file 00515 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info, 00516 mc_molecule_info_test 00517 TYPE(particle_list_p_type), 00518 DIMENSION(:), POINTER :: particles_old 00519 TYPE(particle_list_type), POINTER :: particles_insert, 00520 particles_remove 00521 00522 ! begin the timing of the subroutine 00523 00524 CALL timeset(routineN,handle) 00525 00526 ! reset the overlap flag 00527 loverlap=.FALSE. 00528 00529 ! nullify some pointers 00530 NULLIFY(particles_old,mol_type,mol_type_test,mc_input_file,mc_bias_file) 00531 NULLIFY(oldsys,atom_names,pmswap_mol,insert_coords,remove_coords) 00532 NULLIFY(eta_insert,eta_remove) 00533 00534 ! grab some stuff from mc_par 00535 CALL get_mc_par(mc_par(1)%mc_par,ionode=ionode,BETA=BETA,& 00536 max_val=max_val, min_val=min_val, exp_max_val=exp_max_val,& 00537 exp_min_val=exp_min_val,nswapmoves=nswapmoves,group=group,source=source,& 00538 lbias=lbias,dat_file=dat_file(1),fft_lib=fft_lib,& 00539 mc_molecule_info=mc_molecule_info,pmswap_mol=pmswap_mol) 00540 CALL get_mc_molecule_info(mc_molecule_info,nchains=nchains,& 00541 nunits=nunits,nunits_tot=nunits_tot,nmol_types=nmol_types,& 00542 atom_names=atom_names,mass=mass,mol_type=mol_type) 00543 00544 print_level = 1 00545 00546 CALL get_mc_par(mc_par(2)%mc_par,dat_file=dat_file(2)) 00547 00548 ! allocate some stuff 00549 ALLOCATE(oldsys(1:2),STAT=istat) 00550 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00551 "oldsys",2) 00552 ALLOCATE(particles_old(1:2),STAT=istat) 00553 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00554 "particles_old",2) 00555 00556 ! get the old coordinates 00557 DO ibox=1,2 00558 CALL force_env_get(force_env(ibox)%force_env,& 00559 subsys=oldsys(ibox)%subsys,error=error) 00560 CALL cp_subsys_get(oldsys(ibox)%subsys, & 00561 particles=particles_old(ibox)%list, & 00562 error=error) 00563 ENDDO 00564 00565 ! choose a direction to swap 00566 IF(ionode) rand=next_random_number(rng_stream,error=error) 00567 CALL mp_bcast(rand,source,group) 00568 00569 IF ( rand .LE. 0.50E0_dp ) THEN 00570 remove_box=1 00571 insert_box=2 00572 ELSE 00573 remove_box=2 00574 insert_box=1 00575 ENDIF 00576 00577 ! now assign the eta values for the insert and remove boxes 00578 CALL get_mc_par(mc_par(remove_box)%mc_par,eta=eta_remove) 00579 CALL get_mc_par(mc_par(insert_box)%mc_par,eta=eta_insert) 00580 00581 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! testing 00582 ! remove_box=2 00583 ! insert_box=1 00584 00585 ! now choose a molecule type at random 00586 IF(ionode) rand=next_random_number(rng_stream,error=error) 00587 CALL mp_bcast(rand,source,group) 00588 DO itype=1,nmol_types 00589 IF(rand .LT. pmswap_mol(itype)) THEN 00590 molecule_type=itype 00591 EXIT 00592 ENDIF 00593 ENDDO 00594 00595 ! record the attempt for the box the particle is to be inserted into 00596 moves(molecule_type,insert_box)%moves%swap%attempts=& 00597 moves(molecule_type,insert_box)%moves%swap%attempts+1 00598 00599 ! now choose a random molecule to remove from the removal box, checking 00600 ! to make sure the box isn't empty 00601 IF(nchains(molecule_type,remove_box) == 0) THEN 00602 loverlap=.TRUE. 00603 moves(molecule_type,insert_box)%moves%empty=& 00604 moves(molecule_type,insert_box)%moves%empty+1 00605 ELSE 00606 00607 IF(ionode) rand=next_random_number(rng_stream,error=error) 00608 CALL mp_bcast(rand,source,group) 00609 imolecule=CEILING(rand*nchains(molecule_type,remove_box)) 00610 ! figure out the atom number this molecule starts on 00611 start_atom_rem=1 00612 DO itype=1,nmol_types 00613 IF(itype == molecule_type) THEN 00614 start_atom_rem=start_atom_rem+(imolecule-1)*nunits(itype) 00615 EXIT 00616 ELSE 00617 start_atom_rem=start_atom_rem+nchains(itype,remove_box)*nunits(itype) 00618 ENDIF 00619 ENDDO 00620 00621 ! check for overlap 00622 start_mol=1 00623 DO jbox=1,remove_box-1 00624 start_mol=start_mol+SUM(nchains(:,jbox)) 00625 ENDDO 00626 end_mol=start_mol+SUM(nchains(:,remove_box))-1 00627 CALL check_for_overlap(force_env(remove_box)%force_env,& 00628 nchains(:,remove_box),nunits,loverlap,mol_type(start_mol:end_mol)) 00629 IF(loverlap) CALL stop_program(routineN,moduleN,__LINE__,& 00630 'CBMC swap move found an overlap in the old remove config') 00631 start_mol=1 00632 DO jbox=1,insert_box-1 00633 start_mol=start_mol+SUM(nchains(:,jbox)) 00634 ENDDO 00635 end_mol=start_mol+SUM(nchains(:,insert_box))-1 00636 CALL check_for_overlap(force_env(insert_box)%force_env,& 00637 nchains(:,insert_box),nunits,loverlap,mol_type(start_mol:end_mol)) 00638 IF(loverlap) CALL stop_program(routineN,moduleN,__LINE__,& 00639 'CBMC swap move found an overlap in the old insert config') 00640 ENDIF 00641 00642 IF(loverlap) THEN 00643 DEALLOCATE(oldsys,STAT=istat) 00644 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00645 __LINE__,"oldsys") 00646 DEALLOCATE(particles_old,STAT=istat) 00647 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00648 __LINE__,"particles_old") 00649 CALL timestop(handle) 00650 RETURN 00651 ENDIF 00652 00653 ! figure out how many atoms will be in each box after the move 00654 ins_atoms=nunits_tot(insert_box)+nunits(molecule_type) 00655 rem_atoms=nunits_tot(remove_box)-nunits(molecule_type) 00656 ! now allocate the arrays that will hold the coordinates and the 00657 ! atom name, for writing to the dat file 00658 IF(rem_atoms == 0) THEN 00659 ALLOCATE(remove_coords(1:3,1:nunits(1)),STAT=istat) 00660 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00661 "remove_coords",nunits(1)*3*dp_size) 00662 ALLOCATE(atom_names_remove(1:nunits(1)),STAT=istat) 00663 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00664 "atom_names_remove",nunits(1)*dp_size) 00665 ELSE 00666 ALLOCATE(remove_coords(1:3,1:rem_atoms),STAT=istat) 00667 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00668 "remove_coords",rem_atoms*3*dp_size) 00669 ALLOCATE(atom_names_remove(1:rem_atoms),STAT=istat) 00670 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00671 "atom_names_remove",rem_atoms*dp_size) 00672 ENDIF 00673 ALLOCATE(insert_coords(1:3,1:ins_atoms),STAT=istat) 00674 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00675 "insert_coords",ins_atoms*3*dp_size) 00676 ALLOCATE(atom_names_insert(1:ins_atoms),STAT=istat) 00677 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00678 "atom_names_insert",MAX(ins_atoms,rem_atoms)*dp_size) 00679 00680 ! grab the cells for later...acceptance and insertion 00681 IF(lbias) THEN 00682 CALL force_env_get(bias_env(insert_box)%force_env,& 00683 cell=cell_insert,error=error) 00684 CALL force_env_get(bias_env(remove_box)%force_env,& 00685 cell=cell_remove,error=error) 00686 ELSE 00687 CALL force_env_get(force_env(insert_box)%force_env,& 00688 cell=cell_insert,error=error) 00689 CALL force_env_get(force_env(remove_box)%force_env,& 00690 cell=cell_remove,error=error) 00691 ENDIF 00692 CALL get_cell(cell_remove,abc=abc_remove) 00693 CALL get_cell(cell_insert,abc=abc_insert) 00694 00695 IF(ionode) THEN 00696 ! choose an insertion point 00697 DO idim=1,3 00698 rand=next_random_number(rng_stream,error=error) 00699 pos_insert(idim)=rand*abc_insert(1) 00700 ENDDO 00701 ENDIF 00702 CALL mp_bcast(pos_insert,source,group) 00703 00704 ! allocate some arrays we'll be using 00705 ALLOCATE(r_insert_mol(1:3,1:nunits(molecule_type)),STAT=istat) 00706 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00707 "r_insert_mol",3*nunits(molecule_type)) 00708 00709 iiatom=1 00710 DO iatom=start_atom_rem,start_atom_rem+nunits(molecule_type)-1 00711 r_insert_mol(1:3,iiatom)=& 00712 particles_old(remove_box)%list%els(iatom)%r(1:3) 00713 iiatom=iiatom+1 00714 ENDDO 00715 00716 ! find the center of mass of the molecule 00717 CALL get_center_of_mass(r_insert_mol(:,:),nunits(molecule_type),& 00718 center_of_mass(:),mass(:,molecule_type)) 00719 00720 ! move the center of mass to the insertion point 00721 displace_molecule(1:3)=pos_insert(1:3)-center_of_mass(1:3) 00722 DO iatom=1,nunits(molecule_type) 00723 r_insert_mol(1:3,iatom)=r_insert_mol(1:3,iatom)+& 00724 displace_molecule(1:3) 00725 ENDDO 00726 00727 ! prepare the insertion coordinates to be written to the .dat file so 00728 ! we can create a new force environment...remember there is still a particle 00729 ! in the box even if nchain=0 00730 IF (SUM(nchains(:,insert_box)) == 0) THEN 00731 DO iatom=1,nunits(molecule_type) 00732 insert_coords(1:3,iatom)=r_insert_mol(1:3,iatom) 00733 atom_names_insert(iatom)=& 00734 particles_old(remove_box)%list%els(start_atom_rem+iatom-1)%atomic_kind%name 00735 ENDDO 00736 start_atom_ins=1 00737 ELSE 00738 ! the problem is I can't just tack the new molecule on to the end, 00739 ! because of reading in the dat_file...the topology stuff will crash 00740 ! if the molecules aren't all grouped together, so I have to insert it 00741 ! at the end of the section of molecules with the same type, then 00742 ! remember the start number for the CBMC stuff 00743 start_atom_ins=1 00744 DO itype=1,nmol_types 00745 start_atom_ins=start_atom_ins+& 00746 nchains(itype,insert_box)*nunits(itype) 00747 IF(itype == molecule_type) EXIT 00748 ENDDO 00749 00750 DO iatom=1,start_atom_ins-1 00751 insert_coords(1:3,iatom)=& 00752 particles_old(insert_box)%list%els(iatom)%r(1:3) 00753 atom_names_insert(iatom)=& 00754 particles_old(insert_box)%list%els(iatom)%atomic_kind%name 00755 ENDDO 00756 iiatom=1 00757 DO iatom=start_atom_ins,start_atom_ins+nunits(molecule_type)-1 00758 insert_coords(1:3,iatom)=r_insert_mol(1:3,iiatom) 00759 atom_names_insert(iatom)=atom_names(iiatom,molecule_type) 00760 iiatom=iiatom+1 00761 ENDDO 00762 DO iatom=start_atom_ins+nunits(molecule_type),ins_atoms 00763 insert_coords(1:3,iatom)=& 00764 particles_old(insert_box)%list%els(iatom-nunits(molecule_type))%r(1:3) 00765 atom_names_insert(iatom)=& 00766 particles_old(insert_box)%list%els(iatom-nunits(molecule_type))%atomic_kind%name 00767 ENDDO 00768 ENDIF 00769 00770 ! fold the coordinates into the box and check for overlaps 00771 start_mol=1 00772 DO jbox=1,insert_box-1 00773 start_mol=start_mol+SUM(nchains(:,jbox)) 00774 ENDDO 00775 end_mol=start_mol+SUM(nchains(:,insert_box))-1 00776 00777 ! make the .dat file 00778 IF (ionode) THEN 00779 00780 nchains(molecule_type,insert_box)=nchains(molecule_type,insert_box)+1 00781 IF(lbias) THEN 00782 CALL get_mc_par(mc_par(insert_box)%mc_par,mc_bias_file=mc_bias_file) 00783 CALL mc_make_dat_file_new(insert_coords(:,:),atom_names_insert(:),ins_atoms,& 00784 abc_insert(:),dat_file(insert_box),nchains(:,insert_box),& 00785 mc_bias_file) 00786 ELSE 00787 CALL get_mc_par(mc_par(insert_box)%mc_par,mc_input_file=mc_input_file) 00788 CALL mc_make_dat_file_new(insert_coords(:,:),atom_names_insert(:),ins_atoms,& 00789 abc_insert(:),dat_file(insert_box),nchains(:,insert_box),& 00790 mc_input_file) 00791 ENDIF 00792 nchains(molecule_type,insert_box)=nchains(molecule_type,insert_box)-1 00793 00794 ENDIF 00795 00796 ! now do the same for the removal box...be careful not to make an empty box 00797 IF (rem_atoms == 0) THEN 00798 DO iatom=1,nunits(molecule_type) 00799 remove_coords(1:3,iatom)=r_insert_mol(1:3,iatom) 00800 atom_names_remove(iatom)=atom_names(iatom,molecule_type) 00801 ENDDO 00802 00803 ! need to adjust nchains, because otherwise if we are removing a molecule type 00804 ! that is not the first molecule, the dat file will have two molecules in it but 00805 ! only the coordinates for one 00806 nchains(molecule_type,remove_box)=nchains(molecule_type,remove_box)-1 00807 IF(ionode) THEN 00808 IF(lbias) THEN 00809 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_bias_file=mc_bias_file) 00810 CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 00811 abc_remove(:),dat_file(remove_box),nchains(:,remove_box),& 00812 mc_bias_file) 00813 ELSE 00814 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_input_file=mc_input_file) 00815 CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 00816 abc_remove(:),dat_file(remove_box),nchains(:,remove_box),& 00817 mc_input_file) 00818 ENDIF 00819 00820 ENDIF 00821 nchains(molecule_type,remove_box)=nchains(molecule_type,remove_box)+1 00822 00823 ELSE 00824 DO iatom=1,start_atom_rem-1 00825 remove_coords(1:3,iatom)=& 00826 particles_old(remove_box)%list%els(iatom)%r(1:3) 00827 atom_names_remove(iatom)=& 00828 particles_old(remove_box)%list%els(iatom)%atomic_kind%name 00829 ENDDO 00830 DO iatom=start_atom_rem+nunits(molecule_type),nunits_tot(remove_box) 00831 remove_coords(1:3,iatom-nunits(molecule_type))=& 00832 particles_old(remove_box)%list%els(iatom)%r(1:3) 00833 atom_names_remove(iatom-nunits(molecule_type))=& 00834 particles_old(remove_box)%list%els(iatom)%atomic_kind%name 00835 ENDDO 00836 00837 ! make the .dat file 00838 IF (ionode) THEN 00839 nchains(molecule_type,remove_box)=nchains(molecule_type,remove_box)-1 00840 IF(lbias) THEN 00841 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_bias_file=mc_bias_file) 00842 CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 00843 abc_remove(:),dat_file(remove_box),nchains(:,remove_box),& 00844 mc_bias_file) 00845 ELSE 00846 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_input_file=mc_input_file) 00847 CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 00848 abc_remove(:),dat_file(remove_box),nchains(:,remove_box),& 00849 mc_input_file) 00850 ENDIF 00851 nchains(molecule_type,remove_box)=nchains(molecule_type,remove_box)+1 00852 00853 ENDIF 00854 ENDIF 00855 00856 ! deallocate r_insert_mol 00857 DEALLOCATE(r_insert_mol,STAT=istat) 00858 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00859 "r_insert_mol") 00860 00861 ! now let's create the two new environments with the different number 00862 ! of molecules 00863 ALLOCATE(test_env(1:2),STAT=istat) 00864 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00865 "test_env",3*ins_atoms*dp_size) 00866 CALL mc_create_force_env(test_env(insert_box)%force_env, para_env, & 00867 dat_file(insert_box),error=error) 00868 CALL mc_create_force_env(test_env(remove_box)%force_env, para_env, & 00869 dat_file(remove_box),error=error) 00870 00871 ! allocate an array we'll need 00872 ALLOCATE(r_cbmc(1:3,1:ins_atoms),STAT=istat) 00873 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00874 "r_cbmc",3*ins_atoms*dp_size) 00875 ALLOCATE(cbmc_energies(1:nswapmoves,1:2),STAT=istat) 00876 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00877 "cbmc_energies",nswapmoves*dp_size) 00878 00879 loverlap_ins=.FALSE. 00880 loverlap_rem=.FALSE. 00881 00882 ! compute the new molecule information...we need this for the CBMC part 00883 IF(rem_atoms == 0) THEN 00884 CALL mc_determine_molecule_info(test_env,mc_molecule_info_test,error,& 00885 box_number=remove_box) 00886 ELSE 00887 CALL mc_determine_molecule_info(test_env,mc_molecule_info_test,error) 00888 ENDIF 00889 CALL get_mc_molecule_info(mc_molecule_info_test,nchains=nchains_test,& 00890 mol_type=mol_type_test) 00891 00892 ! figure out the position of the molecule we're inserting, and the 00893 ! Rosenbluth weight 00894 start_mol=1 00895 DO jbox=1,insert_box-1 00896 start_mol=start_mol+SUM(nchains_test(:,jbox)) 00897 ENDDO 00898 end_mol=start_mol+SUM(nchains_test(:,insert_box))-1 00899 00900 IF(lbias) THEN 00901 CALL generate_cbmc_swap_config(test_env(insert_box)%force_env,& 00902 BETA,max_val, min_val, exp_max_val,& 00903 exp_min_val,nswapmoves,weight_new,start_atom_ins,ins_atoms,nunits(:),& 00904 nunits(molecule_type),mass(:,molecule_type),loverlap_ins,bias_energy_new(insert_box),& 00905 bias_energy_old(insert_box),ionode,.FALSE.,mol_type_test(start_mol:end_mol),& 00906 nchains_test(:,insert_box),source,group,rng_stream,error) 00907 00908 ! the energy that comes out of the above routine is the difference...we want 00909 ! the real energy for the acceptance rule...we don't do this for the 00910 ! lbias=.false. case because it doesn't appear in the acceptance rule, and 00911 ! we compensate in case of acceptance 00912 bias_energy_new(insert_box)=bias_energy_new(insert_box)+& 00913 bias_energy_old(insert_box) 00914 ELSE 00915 CALL generate_cbmc_swap_config(test_env(insert_box)%force_env,& 00916 BETA,max_val, min_val, exp_max_val,& 00917 exp_min_val,nswapmoves,weight_new,start_atom_ins,ins_atoms,nunits(:),& 00918 nunits(molecule_type),mass(:,molecule_type),loverlap_ins,new_energy(insert_box),& 00919 old_energy(insert_box),ionode,.FALSE.,mol_type_test(start_mol:end_mol),& 00920 nchains_test(:,insert_box),source,group,rng_stream,error) 00921 ENDIF 00922 00923 CALL force_env_get(test_env(insert_box)%force_env,& 00924 subsys=insert_sys,error=error) 00925 CALL cp_subsys_get(insert_sys, & 00926 particles=particles_insert,error=error) 00927 00928 DO iatom=1,ins_atoms 00929 r_cbmc(1:3,iatom)=particles_insert%els(iatom)%r(1:3) 00930 ENDDO 00931 00932 ! make sure there is no overlap 00933 00934 IF(loverlap_ins .OR. loverlap_rem) THEN 00935 ! deallocate some stuff 00936 CALL mc_molecule_info_destroy(mc_molecule_info_test) 00937 CALL section_vals_release(test_env(insert_box)%force_env%root_section,error=error) 00938 CALL section_vals_release(test_env(remove_box)%force_env%root_section,error=error) 00939 CALL force_env_release(test_env(insert_box)%force_env,error=error) 00940 CALL force_env_release(test_env(remove_box)%force_env,error=error) 00941 DEALLOCATE(insert_coords,STAT=istat) 00942 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00943 __LINE__,"insert_coords") 00944 DEALLOCATE(remove_coords,STAT=istat) 00945 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00946 __LINE__,"remove_coords") 00947 DEALLOCATE(r_cbmc,STAT=istat) 00948 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00949 __LINE__,"r_cbmc") 00950 DEALLOCATE(cbmc_energies,STAT=istat) 00951 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00952 __LINE__,"cbmc_energies") 00953 DEALLOCATE(oldsys,STAT=istat) 00954 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00955 __LINE__,"oldsys") 00956 DEALLOCATE(particles_old,STAT=istat) 00957 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00958 __LINE__,"particles_old") 00959 DEALLOCATE(test_env,STAT=istat) 00960 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 00961 __LINE__,"test_env") 00962 CALL timestop(handle) 00963 RETURN 00964 ENDIF 00965 00966 ! broadcast the choosen coordiantes to all processors 00967 00968 CALL force_env_get(test_env(insert_box)%force_env,& 00969 subsys=insert_sys,error=error) 00970 CALL cp_subsys_get(insert_sys, & 00971 particles=particles_insert,error=error) 00972 00973 DO iatom=1,ins_atoms 00974 particles_insert%els(iatom)%r(1:3)= & 00975 r_cbmc(1:3,iatom) 00976 ENDDO 00977 00978 ! if we made it this far, we have no overlaps 00979 moves(molecule_type,insert_box)%moves%grown=& 00980 moves(molecule_type,insert_box)%moves%grown+1 00981 00982 ! if we're biasing, we need to make environments with the non-biasing 00983 ! potentials, and calculate the energies 00984 IF(lbias) THEN 00985 00986 ALLOCATE(test_env_bias(1:2),STAT=istat) 00987 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00988 "test_env_bias",2) 00989 00990 ! first, the environment to which we added a molecule 00991 CALL get_mc_par(mc_par(insert_box)%mc_par,mc_input_file=mc_input_file) 00992 IF(ionode) CALL mc_make_dat_file_new(r_cbmc(:,:),atom_names_insert(:),ins_atoms,& 00993 abc_insert(:),dat_file(insert_box),nchains_test(:,insert_box),& 00994 mc_input_file) 00995 test_env_bias(insert_box)%force_env => test_env(insert_box)%force_env 00996 NULLIFY(test_env(insert_box)%force_env) 00997 CALL mc_create_force_env(test_env(insert_box)%force_env, para_env, & 00998 dat_file(insert_box),error=error) 00999 ! CALL section_vals_release(insert_env%root_section) 01000 01001 CALL force_env_calc_energy_force(test_env(insert_box)%force_env,& 01002 calc_force=.FALSE.,error=error) 01003 CALL force_env_get(test_env(insert_box)%force_env,& 01004 potential_energy=new_energy(insert_box),error=error) 01005 01006 ! now the environment that has one less molecule 01007 IF (SUM(nchains_test(:,remove_box)) == 0) THEN 01008 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_input_file=mc_input_file) 01009 IF(ionode) CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 01010 abc_remove(:),dat_file(remove_box),nchains_test(:,remove_box),& 01011 mc_input_file) 01012 test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env 01013 NULLIFY(test_env(remove_box)%force_env) 01014 CALL mc_create_force_env(test_env(remove_box)%force_env, para_env, & 01015 dat_file(remove_box),error=error) 01016 new_energy(remove_box)=0.0E0_dp 01017 bias_energy_new(remove_box)=0.0E0_dp 01018 ELSE 01019 CALL get_mc_par(mc_par(remove_box)%mc_par,mc_input_file=mc_input_file) 01020 IF(ionode) CALL mc_make_dat_file_new(remove_coords(:,:),atom_names_remove(:),rem_atoms,& 01021 abc_remove(:),dat_file(remove_box),nchains_test(:,remove_box),& 01022 mc_input_file) 01023 test_env_bias(remove_box)%force_env => test_env(remove_box)%force_env 01024 NULLIFY(test_env(remove_box)%force_env) 01025 CALL mc_create_force_env(test_env(remove_box)%force_env, para_env, & 01026 dat_file(remove_box),error=error) 01027 CALL force_env_calc_energy_force(test_env(remove_box)%force_env,& 01028 calc_force=.FALSE.,error=error) 01029 CALL force_env_get(test_env(remove_box)%force_env,& 01030 potential_energy=new_energy(remove_box),error=error) 01031 CALL force_env_calc_energy_force(test_env_bias(remove_box)%force_env,& 01032 calc_force=.FALSE.,error=error) 01033 CALL force_env_get(test_env_bias(remove_box)%force_env,& 01034 potential_energy=bias_energy_new(remove_box),error=error) 01035 ENDIF 01036 ELSE 01037 IF (SUM(nchains_test(:,remove_box)) == 0) THEN 01038 new_energy(remove_box)=0.0E0_dp 01039 ELSE 01040 CALL force_env_calc_energy_force(test_env(remove_box)%force_env,& 01041 calc_force=.FALSE.,error=error) 01042 CALL force_env_get(test_env(remove_box)%force_env,& 01043 potential_energy=new_energy(remove_box),error=error) 01044 ENDIF 01045 ENDIF 01046 01047 ! now we need to figure out the rosenbluth weight for the old configuration... 01048 ! we wait until now to do that because we need the energy of the box that 01049 ! has had a molecule removed...notice we use the environment that has not 01050 ! had a molecule removed for the CBMC configurations, and therefore nchains 01051 ! and mol_type instead of nchains_test and mol_type_test 01052 start_mol=1 01053 DO jbox=1,remove_box-1 01054 start_mol=start_mol+SUM(nchains(:,jbox)) 01055 ENDDO 01056 end_mol=start_mol+SUM(nchains(:,remove_box))-1 01057 IF(lbias) THEN 01058 CALL generate_cbmc_swap_config(bias_env(remove_box)%force_env,& 01059 BETA,max_val, min_val, exp_max_val,& 01060 exp_min_val,nswapmoves,weight_old,start_atom_rem,nunits_tot(remove_box),& 01061 nunits(:),nunits(molecule_type),mass(:,molecule_type),loverlap_rem,rdum,& 01062 bias_energy_new(remove_box),ionode,.TRUE.,mol_type(start_mol:end_mol),& 01063 nchains(:,remove_box),source,group,rng_stream,error) 01064 ELSE 01065 CALL generate_cbmc_swap_config(force_env(remove_box)%force_env,& 01066 BETA,max_val, min_val, exp_max_val,& 01067 exp_min_val,nswapmoves,weight_old,start_atom_rem,nunits_tot(remove_box),& 01068 nunits(:),nunits(molecule_type),mass(:,molecule_type),loverlap_rem,rdum,& 01069 new_energy(remove_box),ionode,.TRUE.,mol_type(start_mol:end_mol),& 01070 nchains(:,remove_box),source,group,rng_stream,error) 01071 ENDIF 01072 01073 ! figure out the prefactor to the boltzmann weight in the acceptance 01074 ! rule, based on numbers of particles and volumes 01075 01076 prefactor=REAL(nchains(molecule_type,remove_box),dp)/ 01077 REAL(nchains(molecule_type,insert_box)+1,dp)* 01078 abc_insert(1)**3/abc_remove(1)**3 01079 01080 IF(lbias) THEN 01081 01082 del_quickstep_energy=(-BETA)*(new_energy(insert_box)-& 01083 old_energy(insert_box)+new_energy(remove_box)-& 01084 old_energy(remove_box)-(bias_energy_new(insert_box)+& 01085 bias_energy_new(remove_box)-bias_energy_old(insert_box)& 01086 -bias_energy_old(remove_box))) 01087 01088 IF (del_quickstep_energy .GT. exp_max_val) THEN 01089 del_quickstep_energy=max_val 01090 ELSEIF(del_quickstep_energy .LT. exp_min_val) THEN 01091 del_quickstep_energy=min_val 01092 ELSE 01093 del_quickstep_energy=EXP(del_quickstep_energy) 01094 ENDIF 01095 w=prefactor*del_quickstep_energy*weight_new/weight_old & 01096 *EXP(BETA*(eta_remove(molecule_type)-eta_insert(molecule_type))) 01097 01098 ELSE 01099 w=prefactor*weight_new/weight_old & 01100 *EXP(BETA*(eta_remove(molecule_type)-eta_insert(molecule_type))) 01101 01102 ENDIF 01103 01104 ! check if the move is accepted 01105 IF(w .GE. 1.0E0_dp) THEN 01106 rand=0.0E0_dp 01107 ELSE 01108 IF(ionode) rand=next_random_number(rng_stream,error=error) 01109 CALL mp_bcast(rand,source,group) 01110 ENDIF 01111 01112 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01113 IF ( rand .LT. w ) THEN 01114 01115 ! accept the move 01116 01117 ! accept the move 01118 moves(molecule_type,insert_box)%moves%swap%successes=& 01119 moves(molecule_type,insert_box)%moves%swap%successes+1 01120 01121 ! we need to compensate for the fact that we take the difference in 01122 ! generate_cbmc_config to keep the exponetials small 01123 IF(.NOT. lbias) THEN 01124 new_energy(insert_box)=new_energy(insert_box)+& 01125 old_energy(insert_box) 01126 ENDIF 01127 01128 DO ibox=1,2 01129 ! update energies 01130 energy_check(ibox)=energy_check(ibox)+(new_energy(ibox)-& 01131 old_energy(ibox)) 01132 old_energy(ibox)=new_energy(ibox) 01133 ! if we're biasing the update the biasing energy 01134 IF (lbias) THEN 01135 last_bias_energy(ibox)=bias_energy_new(ibox) 01136 bias_energy_old(ibox)=bias_energy_new(ibox) 01137 ENDIF 01138 01139 ENDDO 01140 01141 ! change particle numbers...basically destory the old mc_molecule_info and attach 01142 ! the new stuff to the mc_pars 01143 ! figure out the molecule information for the new environments 01144 CALL mc_molecule_info_destroy(mc_molecule_info) 01145 CALL set_mc_par(mc_par(insert_box)%mc_par,mc_molecule_info=mc_molecule_info_test) 01146 CALL set_mc_par(mc_par(remove_box)%mc_par,mc_molecule_info=mc_molecule_info_test) 01147 01148 ! update coordinates 01149 CALL force_env_get(test_env(insert_box)%force_env,& 01150 subsys=insert_sys,error=error) 01151 CALL cp_subsys_get(insert_sys, & 01152 particles=particles_insert,error=error) 01153 DO ipart=1,ins_atoms 01154 r_old(1:3,ipart,insert_box)=particles_insert%els(ipart)%r(1:3) 01155 ENDDO 01156 CALL force_env_get(test_env(remove_box)%force_env,& 01157 subsys=remove_sys,error=error) 01158 CALL cp_subsys_get(remove_sys, & 01159 particles=particles_remove,error=error) 01160 DO ipart=1,rem_atoms 01161 r_old(1:3,ipart,remove_box)=particles_remove%els(ipart)%r(1:3) 01162 ENDDO 01163 01164 ! Only the root_section of the second subforce_env needs to be released 01165 ! The other is pointing to the main force_env and must not be released 01166 01167 ! insertion box 01168 IF (insert_box/=1) THEN 01169 CALL section_vals_release(test_env(insert_box)%force_env%root_section,error=error) 01170 END IF 01171 CALL section_vals_release(force_env(insert_box)%force_env%root_section,error=error) 01172 CALL force_env_release(force_env(insert_box)%force_env,error=error) 01173 force_env(insert_box)%force_env => test_env(insert_box)%force_env 01174 01175 ! removal box 01176 IF (remove_box/=1) THEN 01177 CALL section_vals_release(test_env(remove_box)%force_env%root_section,error=error) 01178 END IF 01179 CALL section_vals_release(force_env(remove_box)%force_env%root_section,error=error) 01180 CALL force_env_release(force_env(remove_box)%force_env,error=error) 01181 force_env(remove_box)%force_env => test_env(remove_box)%force_env 01182 01183 ! if we're biasing, update the bias_env 01184 IF(lbias) THEN 01185 CALL section_vals_release(bias_env(insert_box)%force_env%root_section,error=error) 01186 CALL section_vals_release(bias_env(remove_box)%force_env%root_section,error=error) 01187 01188 CALL force_env_release(bias_env(insert_box)%force_env,error=error) 01189 bias_env(insert_box)%force_env => test_env_bias(insert_box)%force_env 01190 CALL force_env_release(bias_env(remove_box)%force_env,error=error) 01191 bias_env(remove_box)%force_env => test_env_bias(remove_box)%force_env 01192 DEALLOCATE(test_env_bias,STAT=istat) 01193 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01194 __LINE__,"test_env_bias") 01195 ENDIF 01196 01197 ELSE 01198 01199 ! reject the move 01200 CALL mc_molecule_info_destroy(mc_molecule_info_test) 01201 CALL section_vals_release(test_env(insert_box)%force_env%root_section,error=error) 01202 CALL section_vals_release(test_env(remove_box)%force_env%root_section,error=error) 01203 CALL force_env_release(test_env(insert_box)%force_env,error=error) 01204 CALL force_env_release(test_env(remove_box)%force_env,error=error) 01205 IF(lbias) THEN 01206 CALL section_vals_release(test_env_bias(insert_box)%force_env%root_section,error=error) 01207 CALL section_vals_release(test_env_bias(remove_box)%force_env%root_section,error=error) 01208 CALL force_env_release( test_env_bias(insert_box)%force_env,error=error) 01209 CALL force_env_release( test_env_bias(remove_box)%force_env,error=error) 01210 DEALLOCATE(test_env_bias,STAT=istat) 01211 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01212 __LINE__,"test_env_bias") 01213 ENDIF 01214 ENDIF 01215 01216 ! deallocate some stuff 01217 DEALLOCATE(insert_coords,STAT=istat) 01218 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01219 __LINE__,"insert_coords") 01220 DEALLOCATE(remove_coords,STAT=istat) 01221 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01222 __LINE__,"remove_coords") 01223 DEALLOCATE(test_env,STAT=istat) 01224 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01225 __LINE__,"test_env") 01226 DEALLOCATE(cbmc_energies,STAT=istat) 01227 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01228 __LINE__,"cbmc_energies") 01229 DEALLOCATE(r_cbmc,STAT=istat) 01230 IF (istat /= 0) CALL stop_memory(routineN,moduleN,& 01231 __LINE__,"r_cbmc") 01232 DEALLOCATE(oldsys,STAT=istat) 01233 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01234 "oldsys") 01235 DEALLOCATE(particles_old,STAT=istat) 01236 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01237 "particles_old") 01238 01239 ! end the timing 01240 CALL timestop(handle) 01241 01242 END SUBROUTINE mc_ge_swap_move 01243 01244 ! ***************************************************************************** 01267 SUBROUTINE mc_ge_volume_move ( mc_par,force_env, moves,move_updates,& 01268 nnstep,old_energy,energy_check,r_old,rng_stream,error) 01269 01270 TYPE(mc_simulation_parameters_p_type), 01271 DIMENSION(:), POINTER :: mc_par 01272 TYPE(force_env_p_type), DIMENSION(:), 01273 POINTER :: force_env 01274 TYPE(mc_moves_p_type), DIMENSION(:, :), 01275 POINTER :: moves, move_updates 01276 INTEGER, INTENT(IN) :: nnstep 01277 REAL(KIND=dp), DIMENSION(:), 01278 INTENT(INOUT) :: old_energy, energy_check 01279 REAL(KIND=dp), DIMENSION(:, :, :), 01280 INTENT(INOUT) :: r_old 01281 TYPE(rng_stream_type), POINTER :: rng_stream 01282 TYPE(cp_error_type), INTENT(INOUT) :: error 01283 01284 CHARACTER(len=*), PARAMETER :: routineN = 'mc_ge_volume_move', 01285 routineP = moduleN//':'//routineN 01286 01287 CHARACTER(LEN=200) :: fft_lib 01288 CHARACTER(LEN=40), DIMENSION(1:2) :: dat_file 01289 INTEGER :: cl, end_atom, end_mol, group, handle, iatom, ibox, imolecule, 01290 iside, istat, j, jatom, jbox, max_atoms, molecule_index, molecule_type, 01291 print_level, source, start_atom, start_mol 01292 INTEGER, DIMENSION(:), POINTER :: mol_type, nunits, nunits_tot 01293 INTEGER, DIMENSION(:, :), POINTER :: nchains 01294 LOGICAL :: ionode 01295 LOGICAL, ALLOCATABLE, DIMENSION(:) :: loverlap 01296 LOGICAL, DIMENSION(1:2) :: lempty 01297 REAL(dp), DIMENSION(:, :), POINTER :: mass 01298 REAL(KIND=dp) :: BETA, prefactor, rand, 01299 rmvolume, vol_dis, w 01300 REAL(KIND=dp), ALLOCATABLE, 01301 DIMENSION(:, :, :) :: r 01302 REAL(KIND=dp), DIMENSION(1:2) :: new_energy, volume_new, 01303 volume_old 01304 REAL(KIND=dp), DIMENSION(1:3) :: center_of_mass, 01305 center_of_mass_new, diff 01306 REAL(KIND=dp), DIMENSION(1:3, 1:2) :: abc, new_cell_length, 01307 old_cell_length 01308 REAL(KIND=dp), DIMENSION(1:3, 1:3, 1:2) :: hmat_test 01309 TYPE(cell_p_type), DIMENSION(:), POINTER :: cell, cell_old, cell_test 01310 TYPE(cp_subsys_p_type), DIMENSION(:), 01311 POINTER :: oldsys 01312 TYPE(mc_molecule_info_type), POINTER :: mc_molecule_info 01313 TYPE(particle_list_p_type), 01314 DIMENSION(:), POINTER :: particles_old 01315 01316 ! begin the timing of the subroutine 01317 01318 CALL timeset(routineN,handle) 01319 01320 ! nullify some pointers 01321 NULLIFY(particles_old,cell,oldsys,cell_old,cell_test) 01322 01323 ! get some data from mc_par 01324 CALL get_mc_par(mc_par(1)%mc_par,ionode=ionode,source=source,& 01325 group=group,dat_file=dat_file(1),rmvolume=rmvolume,& 01326 BETA=BETA,cl=cl,fft_lib=fft_lib,& 01327 mc_molecule_info=mc_molecule_info) 01328 CALL get_mc_molecule_info(mc_molecule_info,nunits_tot=nunits_tot,& 01329 mass=mass,nchains=nchains,nunits=nunits,mol_type=mol_type) 01330 01331 print_level = 1 01332 CALL get_mc_par(mc_par(2)%mc_par,dat_file=dat_file(2)) 01333 01334 ! allocate some stuff 01335 max_atoms=MAX(nunits_tot(1),nunits_tot(2)) 01336 ALLOCATE(r(1:3,max_atoms,1:2),STAT=istat) 01337 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01338 "r",2*3*max_atoms*dp_size) 01339 ALLOCATE(oldsys(1:2),STAT=istat) 01340 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01341 "oldsys",2) 01342 ALLOCATE(particles_old(1:2),STAT=istat) 01343 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01344 "particles_old",2) 01345 ALLOCATE(cell(1:2),STAT=istat) 01346 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01347 "cell",2) 01348 ALLOCATE(cell_test(1:2),STAT=istat) 01349 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01350 "cell_test",2) 01351 ALLOCATE(cell_old(1:2),STAT=istat) 01352 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01353 "cell_old",2) 01354 ALLOCATE(loverlap(1:2),STAT=istat) 01355 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01356 "loverlap",2) 01357 01358 ! check for empty boxes...need to be careful because we can't build 01359 ! a force_env with no particles 01360 DO ibox=1,2 01361 lempty(ibox)=.FALSE. 01362 IF(SUM(nchains(:,ibox))==0) THEN 01363 lempty(ibox)=.TRUE. 01364 ENDIF 01365 ENDDO 01366 01367 ! record the attempt 01368 DO ibox=1,2 01369 moves(1,ibox)%moves%volume%attempts= & 01370 moves(1,ibox)%moves%volume%attempts+1 01371 move_updates(1,ibox)%moves%volume%attempts=& 01372 move_updates(1,ibox)%moves%volume%attempts+1 01373 ENDDO 01374 01375 ! now let's grab the cell length and particle positions 01376 DO ibox=1,2 01377 CALL force_env_get(force_env(ibox)%force_env,& 01378 subsys=oldsys(ibox)%subsys,cell=cell(ibox)%cell,error=error) 01379 CALL get_cell(cell(ibox)%cell,abc=abc(:,ibox)) 01380 NULLIFY(cell_old(ibox)%cell) 01381 CALL cell_create(cell_old(ibox)%cell,error=error) 01382 CALL cell_clone(cell(ibox)%cell,cell_old(ibox)%cell,error=error) 01383 CALL cp_subsys_get(oldsys(ibox)%subsys, & 01384 particles=particles_old(ibox)%list, error=error) 01385 01386 ! find the old cell length 01387 old_cell_length(1:3,ibox)=abc(1:3,ibox) 01388 01389 ENDDO 01390 01391 DO ibox=1,2 01392 01393 ! save the old coordiantes 01394 DO iatom=1,nunits_tot(ibox) 01395 r(1:3,iatom,ibox)=particles_old(ibox)%list%els(iatom)%r(1:3) 01396 ENDDO 01397 01398 ENDDO 01399 01400 ! call a random number to figure out how far we're moving 01401 IF(ionode) rand=next_random_number(rng_stream,error=error) 01402 CALL mp_bcast(rand,source,group) 01403 01404 vol_dis=rmvolume*(rand-0.5E0_dp)*2.0E0_dp 01405 01406 ! add to one box, subtract from the other 01407 IF(old_cell_length(1,1)*old_cell_length(2,1)*& 01408 old_cell_length(3,1)+vol_dis .LE. (3.0E0_dp/angstrom)**3) & 01409 CALL stop_program(routineN,moduleN,__LINE__,& 01410 'GE_volume moves are trying to make box 1 smaller than 3') 01411 IF(old_cell_length(1,2)*old_cell_length(2,2)*& 01412 old_cell_length(3,2)+vol_dis .LE. (3.0E0_dp/angstrom)**3) & 01413 CALL stop_program(routineN,moduleN,__LINE__,& 01414 'GE_volume moves are trying to make box 2 smaller than 3') 01415 01416 DO iside=1,3 01417 new_cell_length(iside,1)=(old_cell_length(1,1)**3+& 01418 vol_dis)**(1.0E0_dp/3.0E0_dp) 01419 new_cell_length(iside,2)=(old_cell_length(1,2)**3-& 01420 vol_dis)**(1.0E0_dp/3.0E0_dp) 01421 ENDDO 01422 01423 ! now we need to make the new cells 01424 DO ibox=1,2 01425 hmat_test(:,:,ibox)=0.0e0_dp 01426 hmat_test(1,1,ibox)=new_cell_length(1,ibox) 01427 hmat_test(2,2,ibox)=new_cell_length(2,ibox) 01428 hmat_test(3,3,ibox)=new_cell_length(3,ibox) 01429 NULLIFY (cell_test(ibox)%cell) 01430 CALL cell_create(cell_test(ibox)%cell,hmat=hmat_test(:,:,ibox),& 01431 periodic=cell(ibox)%cell%perd,error=error) 01432 CALL force_env_set_cell(force_env(ibox)%force_env,& 01433 cell_test(ibox)%cell,error=error) 01434 ENDDO 01435 01436 DO ibox=1,2 01437 01438 ! save the coords 01439 DO iatom=1,nunits_tot(ibox) 01440 r(1:3,iatom,ibox)=particles_old(ibox)%list%els(iatom)%r(1:3) 01441 ENDDO 01442 01443 ! now we need to scale the coordinates of all the molecules by the 01444 ! center of mass 01445 start_atom=1 01446 molecule_index=1 01447 DO jbox=1,ibox-1 01448 IF(jbox == ibox) EXIT 01449 molecule_index=molecule_index+SUM(nchains(:,jbox)) 01450 ENDDO 01451 DO imolecule=1,SUM(nchains(:,ibox)) 01452 molecule_type=mol_type(imolecule+molecule_index-1) 01453 IF(imolecule .NE. 1) THEN 01454 start_atom=start_atom+nunits(mol_type(imolecule+molecule_index-2)) 01455 ENDIF 01456 end_atom=start_atom+nunits(molecule_type)-1 01457 01458 ! now find the center of mass 01459 CALL get_center_of_mass(r(:,start_atom:end_atom,ibox),& 01460 nunits(molecule_type),center_of_mass(:),mass(:,molecule_type)) 01461 01462 ! scale the center of mass and determine the vector that points from the 01463 ! old COM to the new one 01464 center_of_mass_new(1:3)=center_of_mass(1:3)*& 01465 new_cell_length(1:3,ibox)/old_cell_length(1:3,ibox) 01466 DO j=1,3 01467 diff(j)=center_of_mass_new(j)-center_of_mass(j) 01468 ! now change the particle positions 01469 DO jatom=start_atom,end_atom 01470 particles_old(ibox)%list%els(jatom)%r(j)=& 01471 particles_old(ibox)%list%els(jatom)%r(j)+diff(j) 01472 ENDDO 01473 01474 ENDDO 01475 ENDDO 01476 01477 ! check for any overlaps we might have 01478 start_mol=1 01479 DO jbox=1,ibox-1 01480 start_mol=start_mol+SUM(nchains(:,jbox)) 01481 ENDDO 01482 end_mol=start_mol+SUM(nchains(:,ibox))-1 01483 CALL check_for_overlap(force_env(ibox)%force_env,& 01484 nchains(:,ibox),nunits,loverlap(ibox),mol_type(start_mol:end_mol),& 01485 cell_length=new_cell_length(:,ibox)) 01486 01487 ENDDO 01488 01489 ! determine the overall energy difference 01490 01491 DO ibox=1,2 01492 IF(loverlap(ibox)) CYCLE 01493 ! remake the force environment and calculate the energy 01494 IF (lempty(ibox)) THEN 01495 new_energy(ibox)=0.0E0_dp 01496 ELSE 01497 01498 CALL force_env_calc_energy_force(force_env(ibox)%force_env,& 01499 calc_force=.FALSE.,error=error) 01500 CALL force_env_get(force_env(ibox)%force_env,& 01501 potential_energy=new_energy(ibox),error=error) 01502 01503 ENDIF 01504 ENDDO 01505 01506 ! accept or reject the move 01507 DO ibox=1,2 01508 volume_new(ibox)=new_cell_length(1,ibox)*& 01509 new_cell_length(2,ibox)*new_cell_length(3,ibox) 01510 volume_old(ibox)=old_cell_length(1,ibox)*& 01511 old_cell_length(2,ibox)*old_cell_length(3,ibox) 01512 ENDDO 01513 prefactor=(volume_new(1)/volume_old(1))**(SUM(nchains(:,1)))*& 01514 (volume_new(2)/volume_old(2))**(SUM(nchains(:,2))) 01515 01516 IF( loverlap(1) .OR. loverlap(2)) THEN 01517 w=0.0E0_dp 01518 ELSE 01519 w=prefactor*EXP(-BETA*& 01520 (new_energy(1)+new_energy(2)-& 01521 old_energy(1)-old_energy(2))) 01522 01523 ENDIF 01524 01525 IF ( w .GE. 1.0E0_dp ) THEN 01526 w=1.0E0_dp 01527 rand=0.0E0_dp 01528 ELSE 01529 IF(ionode) rand=next_random_number(rng_stream,error=error) 01530 CALL mp_bcast(rand,source,group) 01531 ENDIF 01532 01533 IF (rand .LT. w) THEN 01534 01535 ! write cell length, volume, density, and trial displacement to a file 01536 IF(ionode) THEN 01537 01538 WRITE(cl,*) nnstep,new_cell_length(1,1)*& 01539 angstrom,vol_dis*(angstrom)**3,new_cell_length(1,2)*& 01540 angstrom 01541 WRITE(cl,*) nnstep,new_energy(1),& 01542 old_energy(1),new_energy(2),old_energy(2) 01543 WRITE(cl,*) prefactor,w 01544 ENDIF 01545 01546 DO ibox=1,2 01547 ! accept the move 01548 moves(1,ibox)%moves%volume%successes=& 01549 moves(1,ibox)%moves%volume%successes+1 01550 move_updates(1,ibox)%moves%volume%successes=& 01551 move_updates(1,ibox)%moves%volume%successes+1 01552 01553 ! update energies 01554 energy_check(ibox)=energy_check(ibox)+(new_energy(ibox)-& 01555 old_energy(ibox)) 01556 old_energy(ibox)=new_energy(ibox) 01557 01558 ! and the new "old" coordiantes 01559 DO iatom=1,nunits_tot(ibox) 01560 r_old(1:3,iatom,ibox)=& 01561 particles_old(ibox)%list%els(iatom)%r(1:3) 01562 ENDDO 01563 01564 ENDDO 01565 ELSE 01566 01567 ! reject the move 01568 ! write cell length, volume, density, and trial displacement to a file 01569 IF(ionode) THEN 01570 01571 WRITE(cl,*) nnstep,new_cell_length(1,1)*& 01572 angstrom,vol_dis*(angstrom)**3,new_cell_length(1,2)*& 01573 angstrom 01574 WRITE(cl,*) nnstep,new_energy(1),& 01575 old_energy(1),new_energy(2),old_energy(2) 01576 WRITE(cl,*) prefactor,w 01577 01578 ENDIF 01579 01580 ! reset the cell and particle positions 01581 DO ibox=1,2 01582 CALL force_env_set_cell(force_env(ibox)%force_env,& 01583 cell_old(ibox)%cell,error=error) 01584 DO iatom=1,nunits_tot(ibox) 01585 particles_old(ibox)%list%els(iatom)%r(1:3)=r_old(1:3,iatom,ibox) 01586 ENDDO 01587 ENDDO 01588 01589 ENDIF 01590 01591 ! free up some memory 01592 DO ibox=1,2 01593 CALL cell_release(cell_test(ibox)%cell,error=error) 01594 CALL cell_release(cell_old(ibox)%cell,error=error) 01595 ENDDO 01596 DEALLOCATE(r,STAT=istat) 01597 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01598 "r") 01599 DEALLOCATE(oldsys,STAT=istat) 01600 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01601 "oldsys") 01602 DEALLOCATE(particles_old,STAT=istat) 01603 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01604 "particles_old") 01605 DEALLOCATE(cell,STAT=istat) 01606 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01607 "cell") 01608 DEALLOCATE(cell_old,STAT=istat) 01609 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01610 "cell_old") 01611 DEALLOCATE(cell_test,STAT=istat) 01612 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01613 "cell_test") 01614 DEALLOCATE(loverlap,STAT=istat) 01615 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01616 "loverlap") 01617 01618 ! end the timing 01619 CALL timestop(handle) 01620 01621 END SUBROUTINE mc_ge_volume_move 01622 01623 END MODULE mc_ge_moves 01624
1.7.3