(git:58e3e09)
mc_misc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief contains miscellaneous subroutines used in the Monte Carlo runs,
10 !> mostly I/O stuff
11 !> \author MJM
12 ! **************************************************************************************************
13 MODULE mc_misc
14  USE cp_files, ONLY: close_file,&
15  open_file
16  USE force_env_types, ONLY: use_fist_force,&
18  USE kinds, ONLY: default_string_length,&
19  dp
20  USE mathconstants, ONLY: pi
21  USE mc_types, ONLY: &
22  accattempt, get_mc_input_file, get_mc_molecule_info, get_mc_par, mc_averages_type, &
23  mc_input_file_type, mc_molecule_info_type, mc_moves_p_type, mc_moves_type, mc_simpar_type
24  USE physcon, ONLY: angstrom
25 #include "../../base/base_uses.f90"
26 
27  IMPLICIT NONE
28 
29  PRIVATE
30 
33 
34  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mc_misc'
35 
36 CONTAINS
37 
38 ! **************************************************************************************************
39 !> \brief initializes the structure that holds running averages of MC variables
40 !> \param averages the mc_averages strucutre you want to initialize
41 !>
42 !> Suitable for parallel.
43 !> \author MJM
44 ! **************************************************************************************************
45  SUBROUTINE mc_averages_create(averages)
46 
47  TYPE(mc_averages_type), POINTER :: averages
48 
49  CHARACTER(len=*), PARAMETER :: routinen = 'mc_averages_create'
50 
51  INTEGER :: handle
52 
53 ! begin the timing of the subroutine
54 
55  CALL timeset(routinen, handle)
56 
57 ! allocate all the structures...not sure why, but it won't work otherwise
58  ALLOCATE (averages)
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
63 
64 ! end the timing
65  CALL timestop(handle)
66 
67  END SUBROUTINE mc_averages_create
68 
69 ! **************************************************************************************************
70 !> \brief deallocates the structure that holds running averages of MC variables
71 !> \param averages the mc_averages strucutre you want to release
72 !>
73 !> Suitable for parallel.
74 !> \author MJM
75 ! **************************************************************************************************
76  SUBROUTINE mc_averages_release(averages)
77 
78  TYPE(mc_averages_type), POINTER :: averages
79 
80  CHARACTER(len=*), PARAMETER :: routinen = 'mc_averages_release'
81 
82  INTEGER :: handle
83 
84 ! begin the timing of the subroutine
85 
86  CALL timeset(routinen, handle)
87 
88 ! deallocate
89  DEALLOCATE (averages)
90 
91  NULLIFY (averages)
92 
93 ! end the timing
94  CALL timestop(handle)
95 
96  END SUBROUTINE mc_averages_release
97 
98 ! **************************************************************************************************
99 !> \brief writes a bunch of simulation data to the specified unit
100 !> \param mc_par the mc parameters for the simulation
101 !> \param all_moves the structure that holds data on how many moves are
102 !> accepted/rejected
103 !> \param iw the unit to write to
104 !> \param energy_check the sum of the energy changes of each move
105 !> \param initial_energy the initial unbiased energy of the system
106 !> \param final_energy the final unbiased energy of the system
107 !> \param averages the structure that holds computed average properties for
108 !> the simulation
109 !>
110 !> Only use in serial.
111 !> \author MJM
112 ! **************************************************************************************************
113  SUBROUTINE final_mc_write(mc_par, all_moves, iw, energy_check, initial_energy, &
114  final_energy, averages)
115 
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, &
120  final_energy
121  TYPE(mc_averages_type), POINTER :: averages
122 
123  CHARACTER(len=*), PARAMETER :: routinen = 'final_mc_write'
124 
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
128  LOGICAL :: lbias
129  REAL(dp), DIMENSION(:), POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
130  rmtrans
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
134 
135 ! begin the timing of the subroutine
136 
137  CALL timeset(routinen, handle)
138 
139  NULLIFY (mc_molecule_info, rmbond, rmangle, rmdihedral, rmrot, rmtrans)
140 
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)
145  CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types)
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))"
149 
150 ! write out some data averaged over the whole simulation
151  WRITE (iw, *)
152  WRITE (iw, '(A,A)') '*****************************************************', &
153  '***************************'
154  WRITE (iw, '(A,T66,F15.8)') "Average Energy [Hartrees]:", &
155  averages%ave_energy
156  IF (pmswap .GT. 0.0e0_dp) THEN
157  WRITE (iw, '(A,T66,F15.8)') "Average number of molecules:", &
158  averages%molecules
159  END IF
160  WRITE (iw, '(A,A,T65,F16.6)') "Average Volume ", &
161  "[angstroms**3]:", averages%ave_volume*angstrom**3
162 
163  WRITE (iw, *)
164 
165 ! write out acceptance rates for the moves
166 
167 ! volume moves
168  WRITE (iw, '(A,A)') '-----------------------------------------------------', &
169  '---------------------------'
170  string2 = "Attempted Accepted Percent"
171  string1 = "Volume Moves"
172  string3 = "Maximum volume displacement [angstroms**3]= "
173  rmvolume = rmvolume*angstrom**3
174  CALL final_move_write(all_moves(1)%moves%volume, string1, string2, iw, &
175  displacement=rmvolume, lbias=.false., format_string=format_string, &
176  string3=string3)
177 
178  IF (pmcltrans .GT. 0.0e0_dp) THEN
179 
180 ! Cluster translation moves
181  WRITE (iw, '(A,A)') '-----------------------------------------------------', &
182  '---------------------------'
183  string2 = "Attempted Accepted Percent"
184  string1 = "Cluster Translation Moves"
185  string3 = "Maximum cluster translational displacement [angstroms]= "
186  rmcltrans = rmcltrans*angstrom
187  CALL final_move_write(all_moves(1)%moves%cltrans, string1, string2, iw, &
188  displacement=rmcltrans, lbias=lbias, format_string=format_string, &
189  string3=string3)
190 
191  IF (lbias) THEN
192  WRITE (iw, '(A)') "Biased Move Data for cluster translation"
193  WRITE (iw, '(A,A)') '-------------------------------------------------', &
194  '-------------------------------'
195 ! Cluster bias translation moves
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, &
201  string3=string3)
202  END IF
203 
204  END IF
205 
206 ! Hybrid MC moves (basically short MD runs)
207  string2 = "Attempted Accepted Percent"
208  string1 = "HMC Moves"
209  CALL final_move_write(all_moves(1)%moves%hmc, string1, string2, iw)
210 
211 ! Quickstep moves (a series of moves with one potential, and then corrected for
212 ! by another potential
213  string2 = "Attempted Accepted Percent"
214  string1 = "Quickstep Moves"
215  CALL final_move_write(all_moves(1)%moves%Quickstep, string1, string2, iw)
216 
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  '---------------------------'
223 
224  moves => all_moves(itype)%moves
225 
226 ! AVBMC 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)
236 
237 ! conformation changes
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, &
258  string3=string3)
259 
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, &
265  string3=string3)
266 
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, &
272  string3=string3)
273 
274  WRITE (iw, '(A,A,I5)') "Conformational Moves Rejected Because", &
275  "Box Was Empty: ", moves%empty_conf
276  WRITE (iw, '(A,A)') '-----------------------------------------------', &
277  '--------------------------------'
278  END IF
279 
280 ! translation moves
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, &
286  string3=string3)
287 
288 ! rotation moves
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, &
294  string3=string3)
295 
296 ! swap moves
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, &
302  moves%empty, &
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, &
309  moves%grown, &
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  '--------------------------------'
321  END IF
322 
323 ! now we write out information on the classical moves, if it's
324 ! a classical simulations
325  IF (lbias) THEN
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, &
334  string3=string3)
335 
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, &
340  string3=string3)
341 
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, &
346  string3=string3)
347 
348  ! translation moves
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, &
353  string3=string3)
354 
355 ! rotation moves
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, &
360  string3=string3)
361 
362  END IF
363 
364  END DO
365 
366 ! see if the energies add up properly
367  IF (abs(initial_energy + energy_check - final_energy) .GT. 0.0000001e0_dp) &
368  THEN
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
373  END IF
374  WRITE (iw, '(A,A)') '****************************************************', &
375  '****************************'
376  WRITE (iw, *)
377 
378 ! end the timing
379  CALL timestop(handle)
380 
381  END SUBROUTINE final_mc_write
382 
383 ! **************************************************************************************************
384 !> \brief ...
385 !> \param move_data ...
386 !> \param string1 ...
387 !> \param string2 ...
388 !> \param iw ...
389 !> \param string3 ...
390 !> \param format_string ...
391 !> \param lbias ...
392 !> \param displacement ...
393 ! **************************************************************************************************
394  SUBROUTINE final_move_write(move_data, string1, string2, iw, string3, &
395  format_string, lbias, displacement)
396 
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
404 
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  '---------------------------------'
416  END IF
417  ELSE
418  IF (.NOT. PRESENT(string3) .OR. .NOT. PRESENT(lbias) .OR. &
419  .NOT. PRESENT(displacement)) THEN
420  WRITE (iw, *) 'MISSING FLAGS IN FINAL_MOVE_WRITE'
421  END IF
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  '---------------------------------'
434  END IF
435  END IF
436 
437  END SUBROUTINE final_move_write
438 
439 ! **************************************************************************************************
440 !> \brief writes a new input file that CP2K can read in for when we want
441 !> to change a force env (change molecule number)...this is much simpler
442 !> than the version I had used to have, and also more flexible (in a way).
443 !> It assumes that &CELL comes before &COORDS, and &COORDS comes before
444 !> &TOPOLOGY, and &TOPOLOGY comes before &GLOBAL (which comes before MC).
445 !> It also assumes that you use &MOL_SET in &TOPOLOGY. Still, many fewer
446 !> assumptions than before.
447 !>
448 !> box_length and coordinates should be passed in a.u.
449 !> \param coordinates the coordinates of the atoms in the force_env (a.u.)
450 !> \param atom_names ...
451 !> \param nunits_tot the total number of atoms
452 !> \param box_length the length of all sides of the simulation box (angstrom)
453 !> \param filename the name of the file to write to
454 !> \param nchains ...
455 !> \param mc_input_file ...
456 !> \author MJM
457 !> \note Only use in serial.
458 ! **************************************************************************************************
459  SUBROUTINE mc_make_dat_file_new(coordinates, atom_names, nunits_tot, &
460  box_length, filename, nchains, mc_input_file)
461 
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
469 
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
479 
480 ! open the file
481 
482  CALL open_file(file_name=filename, unit_number=unit, &
483  file_action='WRITE', file_status='REPLACE')
484 
485 ! get all the information from the input_file_type
486  CALL get_mc_input_file(mc_input_file, text=text, cell_row=cell_row, &
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)
495 
496 ! how many molecule types?
497  nmol_types = SIZE(nchains)
498 
499 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
500 !!! !!!
501 !!! WARNING: This code assumes that some sections of the input file are in a certain order. !!!
502 !!! !!!
503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
504 
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))
511  END DO
512  cpassert(mol_set_nmol_row(nmol_types) < force_eval_row_end)
513 
514 ! write the global section, but replace the RUN_TYPE
515  DO iline = global_row_start, run_type_row - 1
516  WRITE (unit, '(A)') trim(text(iline))
517  END DO
518  SELECT CASE (in_use)
519  CASE (use_fist_force)
520  WRITE (unit, '(A)') ' RUN_TYPE ENERGY_FORCE'
521  CASE (use_qs_force)
522  WRITE (unit, '(A)') ' RUN_TYPE ENERGY_FORCE'
523  END SELECT
524  DO iline = run_type_row + 1, global_row_end
525  WRITE (unit, '(A)') trim(text(iline))
526  END DO
527 
528 ! write the motion section without modifications
529  DO iline = motion_row_start, motion_row_end
530  WRITE (unit, '(A)') trim(text(iline))
531  END DO
532 
533 ! write force_eval section up to the cell lengths
534  DO iline = force_eval_row_start, cell_row - 1
535  WRITE (unit, '(A)') trim(text(iline))
536  END DO
537 
538 ! substitute in the current cell lengths
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)
543 
544 ! now write everything until the coordinates
545  DO iline = cell_row + 1, coord_row_start
546  WRITE (unit, '(A)') trim(text(iline))
547  END DO
548 
549 ! we may pass nunits_tot=0, but we should still have coordinates
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
555  END DO
556  ELSE
557  DO iunit = 1, nunits_tot
558  WRITE (unit, '(5X,A,2X,3(F15.10))') &
559  trim(adjustl(atom_names(iunit))), &
560  coordinates(1:3, iunit)*angstrom
561  END DO
562  END IF
563 
564 ! now we need to write the MOL_SET section
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))
569  END DO
570 
571 ! have to print out one molecule, even if it's empty
572  IF (nunits_tot == 0 .AND. itype == 1) THEN
573  WRITE (mol_string, '(I8)') 1
574  ELSE
575  WRITE (mol_string, '(I8)') nchains(itype)
576  END IF
577 
578  line_text = text(mol_set_nmol_row(itype))
579  line_text(mol_set_nmol_column(itype):mol_set_nmol_column(itype) + 9) = &
580  mol_string(1:10)
581  WRITE (unit, '(A)') trim(line_text)
582  start_line = mol_set_nmol_row(itype) + 1
583  END DO
584 
585 ! write remainder of force_eval section
586  DO iline = start_line, force_eval_row_end
587  WRITE (unit, '(A)') trim(text(iline))
588  END DO
589 
590 ! close the file
591  CALL close_file(unit_number=unit)
592 
593  END SUBROUTINE mc_make_dat_file_new
594 END MODULE mc_misc
595 
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: dumpdcd.F:1008
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
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.
Definition: cp_files.F:308
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.
Definition: cp_files.F:119
Interface for the force calculations.
integer, parameter, public use_qs_force
integer, parameter, public use_fist_force
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
contains miscellaneous subroutines used in the Monte Carlo runs, mostly I/O stuff
Definition: mc_misc.F:13
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...
Definition: mc_misc.F:461
subroutine, public mc_averages_release(averages)
deallocates the structure that holds running averages of MC variables
Definition: mc_misc.F:77
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
Definition: mc_misc.F:115
subroutine, public mc_averages_create(averages)
initializes the structure that holds running averages of MC variables
Definition: mc_misc.F:46
holds all the structure types needed for Monte Carlo, except the mc_environment_type
Definition: mc_types.F:15
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)
...
Definition: mc_types.F:554
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
Definition: mc_types.F:279
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: mc_types.F:405
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144