23 mc_input_file_type, mc_molecule_info_type, mc_moves_p_type, mc_moves_type, mc_simpar_type
25 #include "../../base/base_uses.f90"
34 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mc_misc'
47 TYPE(mc_averages_type),
POINTER :: averages
49 CHARACTER(len=*),
PARAMETER :: routinen =
'mc_averages_create'
55 CALL timeset(routinen, handle)
59 averages%ave_energy = 0.0e0_dp
60 averages%ave_energy_squared = 0.0e0_dp
61 averages%ave_volume = 0.0e0_dp
62 averages%molecules = 0.0e0_dp
78 TYPE(mc_averages_type),
POINTER :: averages
80 CHARACTER(len=*),
PARAMETER :: routinen =
'mc_averages_release'
86 CALL timeset(routinen, handle)
114 final_energy, averages)
116 TYPE(mc_simpar_type),
POINTER :: mc_par
117 TYPE(mc_moves_p_type),
DIMENSION(:),
POINTER :: all_moves
118 INTEGER,
INTENT(IN) :: iw
119 REAL(kind=
dp),
INTENT(IN) :: energy_check, initial_energy, &
121 TYPE(mc_averages_type),
POINTER :: averages
123 CHARACTER(len=*),
PARAMETER :: routinen =
'final_mc_write'
125 CHARACTER(LEN=5) :: molecule_string, tab_string
126 CHARACTER(LEN=default_string_length) :: format_string, string1, string2, string3
127 INTEGER :: handle, itype, nmol_types
129 REAL(
dp),
DIMENSION(:),
POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
131 REAL(kind=
dp) :: pmcltrans, pmswap, rmcltrans, rmvolume
132 TYPE(mc_molecule_info_type),
POINTER :: mc_molecule_info
133 TYPE(mc_moves_type),
POINTER :: moves
137 CALL timeset(routinen, handle)
139 NULLIFY (mc_molecule_info, rmbond, rmangle, rmdihedral, rmrot, rmtrans)
141 CALL get_mc_par(mc_par, pmswap=pmswap, rmvolume=rmvolume, &
142 lbias=lbias, rmbond=rmbond, rmangle=rmangle, rmdihedral=rmdihedral, &
143 rmtrans=rmtrans, rmcltrans=rmcltrans, pmcltrans=pmcltrans, rmrot=rmrot, &
144 mc_molecule_info=mc_molecule_info)
146 WRITE (molecule_string,
'(I2)') nmol_types
147 WRITE (tab_string,
'(I4)') 81 - 11*nmol_types
148 format_string =
"(A,T"//trim(adjustl(tab_string))//
","//trim(adjustl(molecule_string))//
"(2X,F9.6))"
152 WRITE (iw,
'(A,A)')
'*****************************************************', &
153 '***************************'
154 WRITE (iw,
'(A,T66,F15.8)')
"Average Energy [Hartrees]:", &
156 IF (pmswap .GT. 0.0e0_dp)
THEN
157 WRITE (iw,
'(A,T66,F15.8)')
"Average number of molecules:", &
160 WRITE (iw,
'(A,A,T65,F16.6)')
"Average Volume ", &
161 "[angstroms**3]:", averages%ave_volume*
angstrom**3
168 WRITE (iw,
'(A,A)')
'-----------------------------------------------------', &
169 '---------------------------'
170 string2 =
"Attempted Accepted Percent"
171 string1 =
"Volume Moves"
172 string3 =
"Maximum volume displacement [angstroms**3]= "
174 CALL final_move_write(all_moves(1)%moves%volume, string1, string2, iw, &
175 displacement=rmvolume, lbias=.false., format_string=format_string, &
178 IF (pmcltrans .GT. 0.0e0_dp)
THEN
181 WRITE (iw,
'(A,A)')
'-----------------------------------------------------', &
182 '---------------------------'
183 string2 =
"Attempted Accepted Percent"
184 string1 =
"Cluster Translation Moves"
185 string3 =
"Maximum cluster translational displacement [angstroms]= "
187 CALL final_move_write(all_moves(1)%moves%cltrans, string1, string2, iw, &
188 displacement=rmcltrans, lbias=lbias, format_string=format_string, &
192 WRITE (iw,
'(A)')
"Biased Move Data for cluster translation"
193 WRITE (iw,
'(A,A)')
'-------------------------------------------------', &
194 '-------------------------------'
196 string2 =
"Attempted Accepted Percent"
197 string1 =
"Cluster Translation Moves"
198 string3 =
"Maximum cluster translational displacement [angstroms]="
199 CALL final_move_write(all_moves(1)%moves%bias_cltrans, string1, string2, iw, &
200 displacement=rmcltrans, lbias=lbias, format_string=format_string, &
207 string2 =
"Attempted Accepted Percent"
208 string1 =
"HMC Moves"
209 CALL final_move_write(all_moves(1)%moves%hmc, string1, string2, iw)
213 string2 =
"Attempted Accepted Percent"
214 string1 =
"Quickstep Moves"
215 CALL final_move_write(all_moves(1)%moves%Quickstep, string1, string2, iw)
217 DO itype = 1, nmol_types
218 WRITE (iw,
'(A,A)')
'-----------------------------------------------------', &
219 '---------------------------'
220 WRITE (iw,
'(A,I5)')
'Move Data for Molecule Type ', itype
221 WRITE (iw,
'(A,A)')
'-----------------------------------------------------', &
222 '---------------------------'
224 moves => all_moves(itype)%moves
227 string2 =
"Attempted Accepted Percent"
228 string1 =
"AVBMC moves from in to in"
229 CALL final_move_write(moves%avbmc_inin, string1, string2, iw)
230 string1 =
"AVBMC moves from in to out"
231 CALL final_move_write(moves%avbmc_inout, string1, string2, iw)
232 string1 =
"AVBMC moves from out to in"
233 CALL final_move_write(moves%avbmc_outin, string1, string2, iw)
234 string1 =
"AVBMC moves from out to out"
235 CALL final_move_write(moves%avbmc_outout, string1, string2, iw)
238 IF (moves%angle%attempts .GT. 0 .OR. &
239 moves%bond%attempts .GT. 0 .OR. &
240 moves%dihedral%attempts .GT. 0)
THEN
241 WRITE (iw,
'(A,T43,A)')
"Conformational Moves", &
242 "Attempted Accepted Percent"
243 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
244 moves%bond%attempts + moves%angle%attempts + &
245 moves%dihedral%attempts, &
246 moves%bond%successes + moves%angle%successes + &
247 moves%dihedral%successes, &
248 REAL(moves%bond%successes + moves%
angle%successes + &
249 moves%dihedral%successes,
dp)/ &
250 REAL(moves%bond%attempts + moves%
angle%attempts + &
251 moves%dihedral%attempts,
dp)*100.0e0_dp
252 string2 =
"Attempted Accepted Percent"
253 string1 =
"Bond Changes"
254 string3 =
"Maximum bond displacement [angstroms]= "
255 rmbond(itype) = rmbond(itype)*
angstrom
256 CALL final_move_write(moves%bond, string1, string2, iw, &
257 displacement=rmbond(itype), lbias=lbias, format_string=format_string, &
260 string1 =
"Angle Changes"
261 string3 =
"Maximum angle displacement [degrees]= "
262 rmangle(itype) = rmangle(itype)/
pi*180.0e0_dp
263 CALL final_move_write(moves%angle, string1, string2, iw, &
264 displacement=rmangle(itype), lbias=lbias, format_string=format_string, &
267 string1 =
"Dihedral Changes"
268 string3 =
"Maximum dihedral displacement [degrees]= "
269 rmdihedral(itype) = rmdihedral(itype)/
pi*180.0e0_dp
270 CALL final_move_write(moves%dihedral, string1, string2, iw, &
271 displacement=rmdihedral(itype), lbias=lbias, format_string=format_string, &
274 WRITE (iw,
'(A,A,I5)')
"Conformational Moves Rejected Because", &
275 "Box Was Empty: ", moves%empty_conf
276 WRITE (iw,
'(A,A)')
'-----------------------------------------------', &
277 '--------------------------------'
281 string1 =
"Translation Moves"
282 string3 =
"Maximum molecular translational displacement [angstroms]= "
283 rmtrans(itype) = rmtrans(itype)*
angstrom
284 CALL final_move_write(moves%trans, string1, string2, iw, &
285 displacement=rmtrans(itype), lbias=lbias, format_string=format_string, &
289 string1 =
"Rotation Moves"
290 string3 =
"Maximum molecular rotational displacement [degrees]= "
291 rmrot(itype) = rmrot(itype)/
pi*180.0e0_dp
292 CALL final_move_write(moves%rot, string1, string2, iw, &
293 displacement=rmrot(itype), lbias=lbias, format_string=format_string, &
297 IF (moves%swap%attempts .GT. 0)
THEN
298 WRITE (iw,
'(A,T43,A)')
"Swap Moves into this box", &
299 "Attempted Empty Percent"
300 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
301 moves%swap%attempts, &
303 REAL(moves%empty,
dp)/ &
304 REAL(moves%swap%attempts,
dp)*100.0e0_dp
305 WRITE (iw,
'(A,T43,A)')
" Growths", &
306 "Attempted Successful Percent"
307 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
308 moves%swap%attempts, &
310 REAL(moves%grown,
dp)/ &
311 REAL(moves%swap%attempts,
dp)*100.0e0_dp
312 WRITE (iw,
'(A,T43,A)')
" Total", &
313 "Attempted Accepted Percent"
314 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
315 moves%swap%attempts, &
316 moves%swap%successes, &
317 REAL(moves%swap%successes,
dp)/ &
318 REAL(moves%swap%attempts,
dp)*100.0e0_dp
319 WRITE (iw,
'(A,A)')
'-----------------------------------------------', &
320 '--------------------------------'
326 WRITE (iw,
'(A)')
"Biased Move Data"
327 WRITE (iw,
'(A,A)')
'-------------------------------------------------', &
328 '-------------------------------'
329 string2 =
"Attempted Accepted Percent"
330 string1 =
"Bond Changes"
331 string3 =
"Maximum bond displacement [angstroms]= "
332 CALL final_move_write(moves%bias_bond, string1, string2, iw, &
333 displacement=rmbond(itype), lbias=lbias, format_string=format_string, &
336 string1 =
"Angle Changes"
337 string3 =
"Maximum angle displacement [degrees]= "
338 CALL final_move_write(moves%bias_angle, string1, string2, iw, &
339 displacement=rmangle(itype), lbias=lbias, format_string=format_string, &
342 string1 =
"Dihedral Changes"
343 string3 =
"Maximum dihedral displacement [degrees]= "
344 CALL final_move_write(moves%bias_dihedral, string1, string2, iw, &
345 displacement=rmdihedral(itype), lbias=lbias, format_string=format_string, &
349 string1 =
"Translation Moves"
350 string3 =
"Maximum molecular translational displacement [angstroms]= "
351 CALL final_move_write(moves%bias_trans, string1, string2, iw, &
352 displacement=rmtrans(itype), lbias=lbias, format_string=format_string, &
356 string1 =
"Rotation Moves"
357 string3 =
"Maximum molecular rotational displacement [degrees]= "
358 CALL final_move_write(moves%bias_rot, string1, string2, iw, &
359 displacement=rmrot(itype), lbias=lbias, format_string=format_string, &
367 IF (abs(initial_energy + energy_check - final_energy) .GT. 0.0000001e0_dp) &
369 WRITE (iw, *)
'!!!!!!! We have an energy problem. !!!!!!!!'
370 WRITE (iw,
'(A,T64,F16.10)')
'Final Energy = ', final_energy
371 WRITE (iw,
'(A,T64,F16.10)')
'Initial Energy + energy_check =', &
372 initial_energy + energy_check
374 WRITE (iw,
'(A,A)')
'****************************************************', &
375 '****************************'
379 CALL timestop(handle)
394 SUBROUTINE final_move_write(move_data, string1, string2, iw, string3, &
395 format_string, lbias, displacement)
397 TYPE(accattempt),
POINTER :: move_data
398 CHARACTER(default_string_length),
INTENT(IN) :: string1, string2
399 INTEGER,
INTENT(IN) :: iw
400 CHARACTER(default_string_length),
INTENT(IN), &
401 OPTIONAL :: string3, format_string
402 LOGICAL,
INTENT(IN),
OPTIONAL :: lbias
403 REAL(
dp),
OPTIONAL :: displacement
405 IF (.NOT.
PRESENT(format_string))
THEN
406 IF (move_data%attempts .GT. 0)
THEN
407 WRITE (iw,
'(A,T43,A)') trim(adjustl(string1)), &
408 trim(adjustl(string2))
409 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
410 move_data%attempts, &
411 move_data%successes, &
412 REAL(move_data%successes,
dp)/ &
413 REAL(move_data%attempts,
dp)*100.0e0_dp
414 WRITE (iw,
'(A,A)')
'-----------------------------------------------', &
415 '---------------------------------'
418 IF (.NOT.
PRESENT(string3) .OR. .NOT.
PRESENT(lbias) .OR. &
419 .NOT.
PRESENT(displacement))
THEN
420 WRITE (iw, *)
'MISSING FLAGS IN FINAL_MOVE_WRITE'
422 IF (move_data%attempts .GT. 0)
THEN
423 WRITE (iw,
'(A,T43,A)') trim(adjustl(string1)), &
424 trim(adjustl(string2))
425 WRITE (iw,
'(T46,I6,9X,I6,7X,F7.3)') &
426 move_data%attempts, &
427 move_data%successes, &
428 REAL(move_data%successes,
dp)/ &
429 REAL(move_data%attempts,
dp)*100.0e0_dp
430 IF (.NOT. lbias)
WRITE (iw,
'(A,T71,F10.5)') &
431 string3, displacement
432 WRITE (iw,
'(A,A)')
'-----------------------------------------------', &
433 '---------------------------------'
437 END SUBROUTINE final_move_write
460 box_length, filename, nchains, mc_input_file)
462 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: coordinates
463 CHARACTER(LEN=*),
DIMENSION(:),
INTENT(IN) :: atom_names
464 INTEGER,
INTENT(IN) :: nunits_tot
465 REAL(kind=
dp),
DIMENSION(1:3),
INTENT(IN) :: box_length
466 CHARACTER(LEN=*),
INTENT(IN) :: filename
467 INTEGER,
DIMENSION(:),
INTENT(IN) :: nchains
468 TYPE(mc_input_file_type),
POINTER :: mc_input_file
470 CHARACTER(60) :: cell_string, mol_string
471 CHARACTER(default_string_length) :: line_text
472 CHARACTER(default_string_length),
DIMENSION(:), &
473 POINTER :: atom_names_empty, text
474 INTEGER :: cell_column, cell_row, coord_row_end, coord_row_start, force_eval_row_end, &
475 force_eval_row_start, global_row_end, global_row_start, iline, in_use, itype, iunit, &
476 motion_row_end, motion_row_start, nmol_types, nunits_empty, run_type_row, start_line, unit
477 INTEGER,
DIMENSION(:),
POINTER :: mol_set_nmol_column, mol_set_nmol_row
478 REAL(
dp),
DIMENSION(:, :),
POINTER :: coordinates_empty
482 CALL open_file(file_name=filename, unit_number=unit, &
483 file_action=
'WRITE', file_status=
'REPLACE')
487 cell_column=cell_column, coord_row_start=coord_row_start, &
488 coord_row_end=coord_row_end, mol_set_nmol_row=mol_set_nmol_row, &
489 mol_set_nmol_column=mol_set_nmol_column, &
490 force_eval_row_start=force_eval_row_start, force_eval_row_end=force_eval_row_end, &
491 global_row_start=global_row_start, global_row_end=global_row_end, &
492 run_type_row=run_type_row, in_use=in_use, atom_names_empty=atom_names_empty, &
493 nunits_empty=nunits_empty, coordinates_empty=coordinates_empty, &
494 motion_row_start=motion_row_start, motion_row_end=motion_row_end)
497 nmol_types =
SIZE(nchains)
505 cpassert(force_eval_row_start < cell_row)
506 cpassert(cell_row < coord_row_start)
507 cpassert(coord_row_start < coord_row_end)
508 cpassert(coord_row_end < mol_set_nmol_row(1))
509 DO itype = 1, nmol_types - 1
510 cpassert(mol_set_nmol_row(itype) < mol_set_nmol_row(itype + 1))
512 cpassert(mol_set_nmol_row(nmol_types) < force_eval_row_end)
515 DO iline = global_row_start, run_type_row - 1
516 WRITE (unit,
'(A)') trim(text(iline))
520 WRITE (unit,
'(A)')
' RUN_TYPE ENERGY_FORCE'
522 WRITE (unit,
'(A)')
' RUN_TYPE ENERGY_FORCE'
524 DO iline = run_type_row + 1, global_row_end
525 WRITE (unit,
'(A)') trim(text(iline))
529 DO iline = motion_row_start, motion_row_end
530 WRITE (unit,
'(A)') trim(text(iline))
534 DO iline = force_eval_row_start, cell_row - 1
535 WRITE (unit,
'(A)') trim(text(iline))
539 WRITE (cell_string,
'(3(F13.8,2X))') box_length(1:3)*
angstrom
540 line_text = text(cell_row)
541 line_text(cell_column:cell_column + 50) = cell_string(1:51)
542 WRITE (unit,
'(A)') trim(line_text)
545 DO iline = cell_row + 1, coord_row_start
546 WRITE (unit,
'(A)') trim(text(iline))
550 IF (nunits_tot == 0)
THEN
551 DO iunit = 1, nunits_empty
552 WRITE (unit,
'(5X,A,2X,3(F15.10))') &
553 trim(adjustl(atom_names_empty(iunit))), &
554 coordinates_empty(1:3, iunit)*
angstrom
557 DO iunit = 1, nunits_tot
558 WRITE (unit,
'(5X,A,2X,3(F15.10))') &
559 trim(adjustl(atom_names(iunit))), &
565 start_line = coord_row_end
566 DO itype = 1, nmol_types
567 DO iline = start_line, mol_set_nmol_row(itype) - 1
568 WRITE (unit,
'(A)') trim(text(iline))
572 IF (nunits_tot == 0 .AND. itype == 1)
THEN
573 WRITE (mol_string,
'(I8)') 1
575 WRITE (mol_string,
'(I8)') nchains(itype)
578 line_text = text(mol_set_nmol_row(itype))
579 line_text(mol_set_nmol_column(itype):mol_set_nmol_column(itype) + 9) = &
581 WRITE (unit,
'(A)') trim(line_text)
582 start_line = mol_set_nmol_row(itype) + 1
586 DO iline = start_line, force_eval_row_end
587 WRITE (unit,
'(A)') trim(text(iline))
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Interface for the force calculations.
integer, parameter, public use_qs_force
integer, parameter, public use_fist_force
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
contains miscellaneous subroutines used in the Monte Carlo runs, mostly I/O stuff
subroutine, public mc_make_dat_file_new(coordinates, atom_names, nunits_tot, box_length, filename, nchains, mc_input_file)
writes a new input file that CP2K can read in for when we want to change a force env (change molecule...
subroutine, public mc_averages_release(averages)
deallocates the structure that holds running averages of MC variables
subroutine, public final_mc_write(mc_par, all_moves, iw, energy_check, initial_energy, final_energy, averages)
writes a bunch of simulation data to the specified unit
subroutine, public mc_averages_create(averages)
initializes the structure that holds running averages of MC variables
holds all the structure types needed for Monte Carlo, except the mc_environment_type
subroutine, public get_mc_molecule_info(mc_molecule_info, nmol_types, nchain_total, nboxes, names, conf_prob, nchains, nunits, mol_type, nunits_tot, in_box, atom_names, mass)
...
subroutine, public get_mc_input_file(mc_input_file, run_type_row, run_type_column, coord_row_start, coord_row_end, cell_row, cell_column, force_eval_row_start, force_eval_row_end, global_row_start, global_row_end, mol_set_nmol_row, mol_set_nmol_column, in_use, text, atom_names_empty, nunits_empty, coordinates_empty, motion_row_end, motion_row_start)
accesses the private elements of the mc_input_file_type
subroutine, public get_mc_par(mc_par, nstep, nvirial, iuptrans, iupcltrans, iupvolume, nmoves, nswapmoves, rm, cl, diff, nstart, source, group, lbias, ionode, lrestart, lstop, rmvolume, rmcltrans, rmbond, rmangle, rmrot, rmtrans, temperature, pressure, rclus, BETA, pmswap, pmvolume, pmtraion, pmtrans, pmcltrans, ensemble, PROGRAM, restart_file_name, molecules_file, moves_file, coords_file, energy_file, displacement_file, cell_file, dat_file, data_file, box2_file, fft_lib, iprint, rcut, ldiscrete, discrete_step, pmavbmc, pbias, avbmc_atom, avbmc_rmin, avbmc_rmax, rmdihedral, input_file, mc_molecule_info, pmswap_mol, pmavbmc_mol, pmtrans_mol, pmrot_mol, pmtraion_mol, mc_input_file, mc_bias_file, pmvol_box, pmclus_box, virial_temps, exp_min_val, exp_max_val, min_val, max_val, eta, pmhmc, pmhmc_box, lhmc, rand2skip)
...
Definition of physical constants:
real(kind=dp), parameter, public angstrom