Line data Source code
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,&
17 : use_qs_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 :
31 : PUBLIC :: final_mc_write, mc_averages_create, mc_averages_release, &
32 : mc_make_dat_file_new
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 26 : 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 26 : CALL timeset(routineN, handle)
56 :
57 : ! allocate all the structures...not sure why, but it won't work otherwise
58 26 : 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 26 : CALL timestop(handle)
66 :
67 26 : 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 26 : 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 26 : CALL timeset(routineN, handle)
87 :
88 : ! deallocate
89 26 : DEALLOCATE (averages)
90 :
91 : NULLIFY (averages)
92 :
93 : ! end the timing
94 26 : CALL timestop(handle)
95 :
96 26 : 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 36 : 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 12 : REAL(dp), DIMENSION(:), POINTER :: rmangle, rmbond, rmdihedral, rmrot, &
130 12 : 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 12 : CALL timeset(routineN, handle)
138 :
139 12 : 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 12 : mc_molecule_info=mc_molecule_info)
145 12 : CALL get_mc_molecule_info(mc_molecule_info, nmol_types=nmol_types)
146 12 : WRITE (molecule_string, '(I2)') nmol_types
147 12 : WRITE (tab_string, '(I4)') 81 - 11*nmol_types
148 12 : 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 12 : WRITE (iw, *)
152 12 : WRITE (iw, '(A,A)') '*****************************************************', &
153 24 : '***************************'
154 12 : WRITE (iw, '(A,T66,F15.8)') "Average Energy [Hartrees]:", &
155 24 : averages%ave_energy
156 12 : IF (pmswap .GT. 0.0E0_dp) THEN
157 2 : WRITE (iw, '(A,T66,F15.8)') "Average number of molecules:", &
158 4 : averages%molecules
159 : END IF
160 12 : WRITE (iw, '(A,A,T65,F16.6)') "Average Volume ", &
161 24 : "[angstroms**3]:", averages%ave_volume*angstrom**3
162 :
163 12 : WRITE (iw, *)
164 :
165 : ! write out acceptance rates for the moves
166 :
167 : ! volume moves
168 12 : WRITE (iw, '(A,A)') '-----------------------------------------------------', &
169 24 : '---------------------------'
170 12 : string2 = "Attempted Accepted Percent"
171 12 : string1 = "Volume Moves"
172 12 : string3 = "Maximum volume displacement [angstroms**3]= "
173 12 : 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 12 : string3=string3)
177 :
178 12 : IF (pmcltrans .GT. 0.0E0_dp) THEN
179 :
180 : ! Cluster translation moves
181 1 : WRITE (iw, '(A,A)') '-----------------------------------------------------', &
182 2 : '---------------------------'
183 1 : string2 = "Attempted Accepted Percent"
184 1 : string1 = "Cluster Translation Moves"
185 1 : string3 = "Maximum cluster translational displacement [angstroms]= "
186 1 : 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 1 : string3=string3)
190 :
191 1 : IF (lbias) THEN
192 0 : WRITE (iw, '(A)') "Biased Move Data for cluster translation"
193 0 : WRITE (iw, '(A,A)') '-------------------------------------------------', &
194 0 : '-------------------------------'
195 : ! Cluster bias translation moves
196 0 : string2 = "Attempted Accepted Percent"
197 0 : string1 = "Cluster Translation Moves"
198 0 : 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 0 : string3=string3)
202 : END IF
203 :
204 : END IF
205 :
206 : ! Hybrid MC moves (basically short MD runs)
207 12 : string2 = "Attempted Accepted Percent"
208 12 : string1 = "HMC Moves"
209 12 : 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 12 : string2 = "Attempted Accepted Percent"
214 12 : string1 = "Quickstep Moves"
215 12 : CALL final_move_write(all_moves(1)%moves%Quickstep, string1, string2, iw)
216 :
217 32 : DO itype = 1, nmol_types
218 20 : WRITE (iw, '(A,A)') '-----------------------------------------------------', &
219 40 : '---------------------------'
220 20 : WRITE (iw, '(A,I5)') 'Move Data for Molecule Type ', itype
221 20 : WRITE (iw, '(A,A)') '-----------------------------------------------------', &
222 40 : '---------------------------'
223 :
224 20 : moves => all_moves(itype)%moves
225 :
226 : ! AVBMC moves
227 20 : string2 = "Attempted Accepted Percent"
228 20 : string1 = "AVBMC moves from in to in"
229 20 : CALL final_move_write(moves%avbmc_inin, string1, string2, iw)
230 20 : string1 = "AVBMC moves from in to out"
231 20 : CALL final_move_write(moves%avbmc_inout, string1, string2, iw)
232 20 : string1 = "AVBMC moves from out to in"
233 20 : CALL final_move_write(moves%avbmc_outin, string1, string2, iw)
234 20 : string1 = "AVBMC moves from out to out"
235 20 : CALL final_move_write(moves%avbmc_outout, string1, string2, iw)
236 :
237 : ! conformation changes
238 : IF (moves%angle%attempts .GT. 0 .OR. &
239 20 : moves%bond%attempts .GT. 0 .OR. &
240 : moves%dihedral%attempts .GT. 0) THEN
241 2 : WRITE (iw, '(A,T43,A)') "Conformational Moves", &
242 4 : "Attempted Accepted Percent"
243 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
244 : moves%bond%attempts + moves%angle%attempts + &
245 2 : moves%dihedral%attempts, &
246 : moves%bond%successes + moves%angle%successes + &
247 2 : 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 4 : moves%dihedral%attempts, dp)*100.0E0_dp
252 2 : string2 = "Attempted Accepted Percent"
253 2 : string1 = "Bond Changes"
254 2 : string3 = "Maximum bond displacement [angstroms]= "
255 2 : 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 2 : string3=string3)
259 :
260 2 : string1 = "Angle Changes"
261 2 : string3 = "Maximum angle displacement [degrees]= "
262 2 : 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 2 : string3=string3)
266 :
267 2 : string1 = "Dihedral Changes"
268 2 : string3 = "Maximum dihedral displacement [degrees]= "
269 2 : 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 2 : string3=string3)
273 :
274 2 : WRITE (iw, '(A,A,I5)') "Conformational Moves Rejected Because", &
275 4 : "Box Was Empty: ", moves%empty_conf
276 2 : WRITE (iw, '(A,A)') '-----------------------------------------------', &
277 4 : '--------------------------------'
278 : END IF
279 :
280 : ! translation moves
281 20 : string1 = "Translation Moves"
282 20 : string3 = "Maximum molecular translational displacement [angstroms]= "
283 20 : 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 20 : string3=string3)
287 :
288 : ! rotation moves
289 20 : string1 = "Rotation Moves"
290 20 : string3 = "Maximum molecular rotational displacement [degrees]= "
291 20 : 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 20 : string3=string3)
295 :
296 : ! swap moves
297 20 : IF (moves%swap%attempts .GT. 0) THEN
298 4 : WRITE (iw, '(A,T43,A)') "Swap Moves into this box", &
299 8 : "Attempted Empty Percent"
300 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
301 4 : moves%swap%attempts, &
302 4 : moves%empty, &
303 : REAL(moves%empty, dp)/ &
304 8 : REAL(moves%swap%attempts, dp)*100.0E0_dp
305 4 : WRITE (iw, '(A,T43,A)') " Growths", &
306 8 : "Attempted Successful Percent"
307 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
308 4 : moves%swap%attempts, &
309 4 : moves%grown, &
310 : REAL(moves%grown, dp)/ &
311 8 : REAL(moves%swap%attempts, dp)*100.0E0_dp
312 4 : WRITE (iw, '(A,T43,A)') " Total", &
313 8 : "Attempted Accepted Percent"
314 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
315 4 : moves%swap%attempts, &
316 4 : moves%swap%successes, &
317 : REAL(moves%swap%successes, dp)/ &
318 8 : REAL(moves%swap%attempts, dp)*100.0E0_dp
319 4 : WRITE (iw, '(A,A)') '-----------------------------------------------', &
320 8 : '--------------------------------'
321 : END IF
322 :
323 : ! now we write out information on the classical moves, if it's
324 : ! a classical simulations
325 32 : IF (lbias) THEN
326 14 : WRITE (iw, '(A)') "Biased Move Data"
327 14 : WRITE (iw, '(A,A)') '-------------------------------------------------', &
328 28 : '-------------------------------'
329 14 : string2 = "Attempted Accepted Percent"
330 14 : string1 = "Bond Changes"
331 14 : 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 14 : string3=string3)
335 :
336 14 : string1 = "Angle Changes"
337 14 : 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 14 : string3=string3)
341 :
342 14 : string1 = "Dihedral Changes"
343 14 : 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 14 : string3=string3)
347 :
348 : ! translation moves
349 14 : string1 = "Translation Moves"
350 14 : 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 14 : string3=string3)
354 :
355 : ! rotation moves
356 14 : string1 = "Rotation Moves"
357 14 : 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 14 : string3=string3)
361 :
362 : END IF
363 :
364 : END DO
365 :
366 : ! see if the energies add up properly
367 12 : IF (ABS(initial_energy + energy_check - final_energy) .GT. 0.0000001E0_dp) &
368 : THEN
369 0 : WRITE (iw, *) '!!!!!!! We have an energy problem. !!!!!!!!'
370 0 : WRITE (iw, '(A,T64,F16.10)') 'Final Energy = ', final_energy
371 0 : WRITE (iw, '(A,T64,F16.10)') 'Initial Energy + energy_check =', &
372 0 : initial_energy + energy_check
373 : END IF
374 12 : WRITE (iw, '(A,A)') '****************************************************', &
375 24 : '****************************'
376 12 : WRITE (iw, *)
377 :
378 : ! end the timing
379 12 : CALL timestop(handle)
380 :
381 12 : 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 233 : 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 233 : IF (.NOT. PRESENT(format_string)) THEN
406 104 : IF (move_data%attempts .GT. 0) THEN
407 7 : WRITE (iw, '(A,T43,A)') TRIM(ADJUSTL(string1)), &
408 14 : TRIM(ADJUSTL(string2))
409 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
410 7 : move_data%attempts, &
411 7 : move_data%successes, &
412 : REAL(move_data%successes, dp)/ &
413 14 : REAL(move_data%attempts, dp)*100.0E0_dp
414 7 : WRITE (iw, '(A,A)') '-----------------------------------------------', &
415 14 : '---------------------------------'
416 : END IF
417 : ELSE
418 129 : IF (.NOT. PRESENT(string3) .OR. .NOT. PRESENT(lbias) .OR. &
419 : .NOT. PRESENT(displacement)) THEN
420 0 : WRITE (iw, *) 'MISSING FLAGS IN FINAL_MOVE_WRITE'
421 : END IF
422 129 : IF (move_data%attempts .GT. 0) THEN
423 54 : WRITE (iw, '(A,T43,A)') TRIM(ADJUSTL(string1)), &
424 108 : TRIM(ADJUSTL(string2))
425 : WRITE (iw, '(T46,I6,9X,I6,7X,F7.3)') &
426 54 : move_data%attempts, &
427 54 : move_data%successes, &
428 : REAL(move_data%successes, dp)/ &
429 108 : REAL(move_data%attempts, dp)*100.0E0_dp
430 54 : IF (.NOT. lbias) WRITE (iw, '(A,T71,F10.5)') &
431 12 : string3, displacement
432 54 : WRITE (iw, '(A,A)') '-----------------------------------------------', &
433 108 : '---------------------------------'
434 : END IF
435 : END IF
436 :
437 233 : 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 102 : SUBROUTINE mc_make_dat_file_new(coordinates, atom_names, nunits_tot, &
460 102 : 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 102 : 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 102 : INTEGER, DIMENSION(:), POINTER :: mol_set_nmol_column, mol_set_nmol_row
478 102 : REAL(dp), DIMENSION(:, :), POINTER :: coordinates_empty
479 :
480 : ! open the file
481 :
482 : CALL open_file(file_name=filename, unit_number=unit, &
483 102 : 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 102 : motion_row_start=motion_row_start, motion_row_end=motion_row_end)
495 :
496 : ! how many molecule types?
497 102 : 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 102 : CPASSERT(force_eval_row_start < cell_row)
506 102 : CPASSERT(cell_row < coord_row_start)
507 102 : CPASSERT(coord_row_start < coord_row_end)
508 102 : CPASSERT(coord_row_end < mol_set_nmol_row(1))
509 203 : DO itype = 1, nmol_types - 1
510 203 : CPASSERT(mol_set_nmol_row(itype) < mol_set_nmol_row(itype + 1))
511 : END DO
512 102 : CPASSERT(mol_set_nmol_row(nmol_types) < force_eval_row_end)
513 :
514 : ! write the global section, but replace the RUN_TYPE
515 408 : DO iline = global_row_start, run_type_row - 1
516 408 : WRITE (unit, '(A)') TRIM(text(iline))
517 : END DO
518 101 : SELECT CASE (in_use)
519 : CASE (use_fist_force)
520 101 : WRITE (unit, '(A)') ' RUN_TYPE ENERGY_FORCE'
521 : CASE (use_qs_force)
522 102 : WRITE (unit, '(A)') ' RUN_TYPE ENERGY_FORCE'
523 : END SELECT
524 204 : DO iline = run_type_row + 1, global_row_end
525 204 : WRITE (unit, '(A)') TRIM(text(iline))
526 : END DO
527 :
528 : ! write the motion section without modifications
529 6028 : DO iline = motion_row_start, motion_row_end
530 6028 : WRITE (unit, '(A)') TRIM(text(iline))
531 : END DO
532 :
533 : ! write force_eval section up to the cell lengths
534 8106 : DO iline = force_eval_row_start, cell_row - 1
535 8106 : WRITE (unit, '(A)') TRIM(text(iline))
536 : END DO
537 :
538 : ! substitute in the current cell lengths
539 408 : WRITE (cell_string, '(3(F13.8,2X))') box_length(1:3)*angstrom
540 102 : line_text = text(cell_row)
541 102 : line_text(cell_column:cell_column + 50) = cell_string(1:51)
542 102 : WRITE (unit, '(A)') TRIM(line_text)
543 :
544 : ! now write everything until the coordinates
545 309 : DO iline = cell_row + 1, coord_row_start
546 309 : WRITE (unit, '(A)') TRIM(text(iline))
547 : END DO
548 :
549 : ! we may pass nunits_tot=0, but we should still have coordinates
550 102 : IF (nunits_tot == 0) THEN
551 0 : DO iunit = 1, nunits_empty
552 : WRITE (unit, '(5X,A,2X,3(F15.10))') &
553 0 : TRIM(ADJUSTL(atom_names_empty(iunit))), &
554 0 : coordinates_empty(1:3, iunit)*angstrom
555 : END DO
556 : ELSE
557 2638 : DO iunit = 1, nunits_tot
558 : WRITE (unit, '(5X,A,2X,3(F15.10))') &
559 2536 : TRIM(ADJUSTL(atom_names(iunit))), &
560 12782 : 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 305 : DO itype = 1, nmol_types
567 1126 : DO iline = start_line, mol_set_nmol_row(itype) - 1
568 1126 : WRITE (unit, '(A)') TRIM(text(iline))
569 : END DO
570 :
571 : ! have to print out one molecule, even if it's empty
572 203 : IF (nunits_tot == 0 .AND. itype == 1) THEN
573 0 : WRITE (mol_string, '(I8)') 1
574 : ELSE
575 203 : WRITE (mol_string, '(I8)') nchains(itype)
576 : END IF
577 :
578 203 : line_text = text(mol_set_nmol_row(itype))
579 : line_text(mol_set_nmol_column(itype):mol_set_nmol_column(itype) + 9) = &
580 203 : mol_string(1:10)
581 203 : WRITE (unit, '(A)') TRIM(line_text)
582 305 : start_line = mol_set_nmol_row(itype) + 1
583 : END DO
584 :
585 : ! write remainder of force_eval section
586 612 : DO iline = start_line, force_eval_row_end
587 612 : WRITE (unit, '(A)') TRIM(text(iline))
588 : END DO
589 :
590 : ! close the file
591 102 : CALL close_file(unit_number=unit)
592 :
593 102 : END SUBROUTINE MC_MAKE_DAT_FILE_NEW
594 : END MODULE mc_misc
595 :
|
