98 #include "./base/base_uses.f90"
103 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'colvar_methods'
104 REAL(KIND=
dp),
PRIVATE,
PARAMETER :: tolerance_acos = 1.0e-5_dp
122 RECURSIVE SUBROUTINE colvar_read(colvar, icol, colvar_section, para_env)
123 TYPE(colvar_type),
POINTER :: colvar
124 INTEGER,
INTENT(IN) :: icol
125 TYPE(section_vals_type),
POINTER :: colvar_section
126 TYPE(mp_para_env_type),
POINTER :: para_env
128 CHARACTER(len=*),
PARAMETER :: routinen =
'colvar_read'
130 CHARACTER(LEN=3) :: fmid
131 CHARACTER(LEN=7) :: tag, tag_comp, tag_comp1, tag_comp2
132 CHARACTER(LEN=default_path_length) :: path_function
133 CHARACTER(LEN=default_string_length) :: tmpstr, tmpstr2
134 CHARACTER(LEN=default_string_length), &
135 DIMENSION(:),
POINTER :: c_kinds, my_par
136 INTEGER :: handle, i, iatm, icomponent, iend, &
137 ifunc, ii, isize, istart, iw, iw1, j, &
138 k, kk, n_var, n_var_k, ncol, ndim, &
140 INTEGER,
DIMENSION(:),
POINTER :: iatms
141 INTEGER,
DIMENSION(:, :),
POINTER :: p_bounds
142 LOGICAL :: check, use_mixed_energy
143 LOGICAL,
DIMENSION(26) :: my_subsection
144 REAL(
dp),
DIMENSION(:),
POINTER :: s1, wei, weights
145 REAL(
dp),
DIMENSION(:, :),
POINTER :: p_range, s1v
146 REAL(kind=
dp),
DIMENSION(1) :: my_val
147 REAL(kind=
dp),
DIMENSION(:),
POINTER :: g_range, grid_point, grid_sp, my_vals, &
149 TYPE(cp_logger_type),
POINTER :: logger
150 TYPE(enumeration_type),
POINTER :: enum
151 TYPE(keyword_type),
POINTER :: keyword
152 TYPE(section_type),
POINTER :: section
153 TYPE(section_vals_type),
POINTER :: acid_hyd_dist_section, acid_hyd_shell_section, &
154 angle_section, colvar_subsection, combine_section, coordination_section, dfunct_section, &
155 distance_from_path_section, distance_section, frame_section, gyration_section, &
156 hbp_section, hydronium_dist_section, hydronium_shell_section, mindist_section, &
157 path_section, plane_dist_section, plane_plane_angle_section, plane_sections, &
158 point_section, population_section, qparm_section, reaction_path_section, &
159 ring_puckering_section, rmsd_section, rotation_section, torsion_section, u_section, &
160 wc_section, wrk_section
161 TYPE(section_vals_type),
POINTER :: xyz_diag_section, xyz_outerdiag_section
163 CALL timeset(routinen, handle)
164 NULLIFY (logger, c_kinds, iatms)
166 my_subsection = .false.
174 plane_plane_angle_section &
181 acid_hyd_shell_section &
184 can_return_null=.true.)
185 distance_from_path_section &
187 i_rep_section=icol, can_return_null=.true.)
189 can_return_null=.true.)
198 ring_puckering_section &
212 IF (
ASSOCIATED(reaction_path_section))
THEN
214 explicit=my_subsection(10))
216 IF (
ASSOCIATED(distance_from_path_section))
THEN
218 explicit=my_subsection(16))
220 IF (
ASSOCIATED(combine_section))
THEN
225 explicit=my_subsection(13))
234 explicit=my_subsection(22))
241 cpassert(count(my_subsection) == 1)
242 cpassert(.NOT.
ASSOCIATED(colvar))
244 IF (my_subsection(1))
THEN
246 wrk_section => distance_section
248 CALL colvar_check_points(colvar, distance_section)
250 colvar%dist_param%i_at = iatms(1)
251 colvar%dist_param%j_at = iatms(2)
254 ELSE IF (my_subsection(2))
THEN
256 wrk_section => angle_section
258 CALL colvar_check_points(colvar, angle_section)
260 colvar%angle_param%i_at_angle = iatms
261 ELSE IF (my_subsection(3))
THEN
263 wrk_section => torsion_section
265 CALL colvar_check_points(colvar, torsion_section)
267 colvar%torsion_param%i_at_tors = iatms
268 colvar%torsion_param%o0 = 0.0_dp
269 ELSE IF (my_subsection(4))
THEN
271 wrk_section => coordination_section
273 CALL colvar_check_points(colvar, coordination_section)
274 NULLIFY (colvar%coord_param%i_at_from, colvar%coord_param%c_kinds_from)
275 NULLIFY (colvar%coord_param%i_at_to, colvar%coord_param%c_kinds_to)
276 NULLIFY (colvar%coord_param%i_at_to_b, colvar%coord_param%c_kinds_to_b)
284 CALL reallocate(colvar%coord_param%i_at_from, 1, ndim +
SIZE(iatms))
285 colvar%coord_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
286 ndim = ndim +
SIZE(iatms)
288 colvar%coord_param%n_atoms_from = ndim
289 colvar%coord_param%use_kinds_from = .false.
296 CALL reallocate(colvar%coord_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
297 colvar%coord_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
298 ndim = ndim +
SIZE(c_kinds)
300 colvar%coord_param%n_atoms_from = 0
301 colvar%coord_param%use_kinds_from = .true.
304 CALL uppercase(colvar%coord_param%c_kinds_from(k))
314 CALL reallocate(colvar%coord_param%i_at_to, 1, ndim +
SIZE(iatms))
315 colvar%coord_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
316 ndim = ndim +
SIZE(iatms)
318 colvar%coord_param%n_atoms_to = ndim
319 colvar%coord_param%use_kinds_to = .false.
326 CALL reallocate(colvar%coord_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
327 colvar%coord_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
328 ndim = ndim +
SIZE(c_kinds)
330 colvar%coord_param%n_atoms_to = 0
331 colvar%coord_param%use_kinds_to = .true.
334 CALL uppercase(colvar%coord_param%c_kinds_to(k))
345 IF (n_var /= 0 .OR. n_var_k /= 0)
THEN
346 colvar%coord_param%do_chain = .true.
351 CALL reallocate(colvar%coord_param%i_at_to_b, 1, ndim +
SIZE(iatms))
352 colvar%coord_param%i_at_to_b(ndim + 1:ndim +
SIZE(iatms)) = iatms
353 ndim = ndim +
SIZE(iatms)
355 colvar%coord_param%n_atoms_to_b = ndim
356 colvar%coord_param%use_kinds_to_b = .false.
360 cpassert(n_var_k > 0)
363 CALL reallocate(colvar%coord_param%c_kinds_to_b, 1, ndim +
SIZE(c_kinds))
364 colvar%coord_param%c_kinds_to_b(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
365 ndim = ndim +
SIZE(c_kinds)
367 colvar%coord_param%n_atoms_to_b = 0
368 colvar%coord_param%use_kinds_to_b = .true.
371 CALL uppercase(colvar%coord_param%c_kinds_to_b(k))
379 colvar%coord_param%do_chain = .false.
380 colvar%coord_param%n_atoms_to_b = 0
381 colvar%coord_param%use_kinds_to_b = .false.
382 NULLIFY (colvar%coord_param%i_at_to_b)
383 NULLIFY (colvar%coord_param%c_kinds_to_b)
384 colvar%coord_param%nncrd_b = 0
385 colvar%coord_param%ndcrd_b = 0
386 colvar%coord_param%r_0_b = 0._dp
389 ELSE IF (my_subsection(5))
THEN
391 wrk_section => plane_dist_section
393 CALL colvar_check_points(colvar, plane_dist_section)
395 cpassert(
SIZE(iatms) == 3)
396 colvar%plane_distance_param%plane = iatms
398 colvar%plane_distance_param%point = iatm
400 ELSE IF (my_subsection(6))
THEN
402 wrk_section => rotation_section
404 CALL colvar_check_points(colvar, rotation_section)
405 CALL section_vals_val_get(rotation_section,
"P1_BOND1", i_val=colvar%rotation_param%i_at1_bond1)
406 CALL section_vals_val_get(rotation_section,
"P2_BOND1", i_val=colvar%rotation_param%i_at2_bond1)
407 CALL section_vals_val_get(rotation_section,
"P1_BOND2", i_val=colvar%rotation_param%i_at1_bond2)
408 CALL section_vals_val_get(rotation_section,
"P2_BOND2", i_val=colvar%rotation_param%i_at2_bond2)
409 ELSE IF (my_subsection(7))
THEN
411 wrk_section => dfunct_section
413 CALL colvar_check_points(colvar, dfunct_section)
415 colvar%dfunct_param%i_at_dfunct = iatms
418 ELSE IF (my_subsection(8))
THEN
420 wrk_section => qparm_section
422 CALL colvar_check_points(colvar, qparm_section)
425 CALL section_vals_val_get(qparm_section,
"INCLUDE_IMAGES", l_val=colvar%qparm_param%include_images)
428 NULLIFY (colvar%qparm_param%i_at_from)
429 NULLIFY (colvar%qparm_param%i_at_to)
434 CALL reallocate(colvar%qparm_param%i_at_from, 1, ndim +
SIZE(iatms))
435 colvar%qparm_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
436 ndim = ndim +
SIZE(iatms)
438 colvar%qparm_param%n_atoms_from = ndim
444 CALL reallocate(colvar%qparm_param%i_at_to, 1, ndim +
SIZE(iatms))
445 colvar%qparm_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
446 ndim = ndim +
SIZE(iatms)
448 colvar%qparm_param%n_atoms_to = ndim
449 ELSE IF (my_subsection(9))
THEN
452 NULLIFY (colvar%hydronium_shell_param%i_oxygens)
453 NULLIFY (colvar%hydronium_shell_param%i_hydrogens)
455 colvar%hydronium_shell_param%n_oxygens, &
456 colvar%hydronium_shell_param%n_hydrogens, &
457 colvar%hydronium_shell_param%i_oxygens, &
458 colvar%hydronium_shell_param%i_hydrogens)
459 ELSE IF (my_subsection(10) .OR. my_subsection(16))
THEN
461 IF (my_subsection(10))
THEN
462 path_section => reaction_path_section
466 ELSE IF (my_subsection(16))
THEN
467 path_section => distance_from_path_section
472 colvar%use_points = .false.
474 CALL section_vals_val_get(path_section,
"DISTANCES_RMSD", l_val=colvar%reaction_path_param%dist_rmsd)
476 IF (colvar%reaction_path_param%dist_rmsd .AND. colvar%reaction_path_param%rmsd)
THEN
477 cpabort(
"CV REACTION PATH: only one between DISTANCES_RMSD and RMSD can be used ")
479 IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd)
THEN
480 NULLIFY (colvar%reaction_path_param%i_rmsd, colvar%reaction_path_param%r_ref)
484 colvar%reaction_path_param%nr_frames = nr_frame
485 CALL read_frames(frame_section, para_env, nr_frame, colvar%reaction_path_param%r_ref, &
486 colvar%reaction_path_param%n_components)
488 IF (colvar%reaction_path_param%subset ==
rmsd_all)
THEN
489 ALLOCATE (colvar%reaction_path_param%i_rmsd(colvar%reaction_path_param%n_components))
490 DO i = 1, colvar%reaction_path_param%n_components
491 colvar%reaction_path_param%i_rmsd(i) = i
493 ELSE IF (colvar%reaction_path_param%subset ==
rmsd_list)
THEN
501 CALL reallocate(colvar%reaction_path_param%i_rmsd, 1, ndim +
SIZE(iatms))
502 colvar%reaction_path_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
503 ndim = ndim +
SIZE(iatms)
505 colvar%reaction_path_param%n_components = ndim
507 cpabort(
"CV REACTION PATH: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
511 CALL section_vals_val_get(path_section,
"ALIGN_FRAMES", l_val=colvar%reaction_path_param%align_frames)
515 ALLOCATE (colvar%reaction_path_param%colvar_p(ncol))
518 NULLIFY (colvar%reaction_path_param%colvar_p(i)%colvar)
519 CALL colvar_read(colvar%reaction_path_param%colvar_p(i)%colvar, i, colvar_subsection, para_env)
522 cpabort(
"CV REACTION PATH: the number of CV to define the path must be >0 ")
524 colvar%reaction_path_param%n_components = ncol
527 CALL section_vals_val_get(path_section,
"STEP_SIZE", r_val=colvar%reaction_path_param%step_size)
528 iend = ceiling(max(range(1), range(2))/colvar%reaction_path_param%step_size)
529 istart = floor(min(range(1), range(2))/colvar%reaction_path_param%step_size)
530 colvar%reaction_path_param%function_bounds(1) = istart
531 colvar%reaction_path_param%function_bounds(2) = iend
532 colvar%reaction_path_param%nr_frames = 2
533 ALLOCATE (colvar%reaction_path_param%f_vals(ncol, istart:iend))
536 check = (ncol ==
SIZE(colvar%reaction_path_param%colvar_p))
541 CALL compress(path_function, full=.true.)
542 CALL parsef(i, trim(path_function), my_par)
544 my_val = real(j, kind=
dp)*colvar%reaction_path_param%step_size
545 colvar%reaction_path_param%f_vals(i, j) =
evalf(i, my_val)
551 "MAP", middle_name=fmid, extension=
".dat", file_status=
"REPLACE")
554 ALLOCATE (grid_sp(ncol))
559 cpassert(ncol ==
SIZE(grid_sp))
560 ALLOCATE (p_range(2, ncol))
561 ALLOCATE (p_bounds(2, ncol))
564 p_range(:, i) = g_range(:)
565 p_bounds(2, i) = ceiling(max(p_range(1, i), p_range(2, i))/grid_sp(i))
566 p_bounds(1, i) = floor(min(p_range(1, i), p_range(2, i))/grid_sp(i))
568 ALLOCATE (s1v(2, istart:iend))
570 ALLOCATE (grid_point(ncol))
572 kk = rec_eval_grid(iw1, ncol, colvar%reaction_path_param%f_vals, v_count, &
573 grid_point, grid_sp, colvar%reaction_path_param%step_size, istart, &
574 iend, s1v, s1, p_bounds, colvar%reaction_path_param%lambda, ifunc=ifunc, &
575 nconf=colvar%reaction_path_param%nr_frames)
578 DEALLOCATE (p_bounds)
581 DEALLOCATE (grid_point)
587 ELSE IF (my_subsection(11))
THEN
590 colvar%use_points = .false.
593 ALLOCATE (colvar%combine_cvs_param%colvar_p(ncol))
596 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
598 WRITE (iw,
'( A )')
' '// &
599 '**********************************************************************'
600 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
601 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| COMBINATION OF THE FOLLOWING COLVARS:'
604 "PRINT%PROGRAM_RUN_INFO")
607 NULLIFY (colvar%combine_cvs_param%colvar_p(i)%colvar)
608 CALL colvar_read(colvar%combine_cvs_param%colvar_p(i)%colvar, i, colvar_subsection, para_env)
612 CALL compress(colvar%combine_cvs_param%function, full=.true.)
615 ALLOCATE (colvar%combine_cvs_param%variables(
SIZE(my_par)))
616 colvar%combine_cvs_param%variables = my_par
618 IF (
SIZE(my_par) /= ncol) &
619 CALL cp_abort(__location__, &
620 "Number of defined COLVAR for COMBINE_COLVAR is different from the "// &
621 "number of variables! It is not possible to define COLVARs in a COMBINE_COLVAR "// &
622 "and avoid their usage in the combininig function!")
624 ALLOCATE (colvar%combine_cvs_param%c_parameters(0))
627 isize =
SIZE(colvar%combine_cvs_param%c_parameters)
629 CALL reallocate(colvar%combine_cvs_param%c_parameters, 1, isize +
SIZE(my_par))
630 colvar%combine_cvs_param%c_parameters(isize + 1:isize +
SIZE(my_par)) = my_par
632 ALLOCATE (colvar%combine_cvs_param%v_parameters(0))
635 isize =
SIZE(colvar%combine_cvs_param%v_parameters)
637 CALL reallocate(colvar%combine_cvs_param%v_parameters, 1, isize +
SIZE(my_vals))
638 colvar%combine_cvs_param%v_parameters(isize + 1:isize +
SIZE(my_vals)) = my_vals
643 ELSE IF (my_subsection(12))
THEN
645 wrk_section => population_section
647 CALL colvar_check_points(colvar, population_section)
649 NULLIFY (colvar%population_param%i_at_from, colvar%population_param%c_kinds_from)
650 NULLIFY (colvar%population_param%i_at_to, colvar%population_param%c_kinds_to)
659 CALL reallocate(colvar%population_param%i_at_from, 1, ndim +
SIZE(iatms))
660 colvar%population_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
661 ndim = ndim +
SIZE(iatms)
663 colvar%population_param%n_atoms_from = ndim
664 colvar%population_param%use_kinds_from = .false.
671 CALL reallocate(colvar%population_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
672 colvar%population_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
673 ndim = ndim +
SIZE(c_kinds)
675 colvar%population_param%n_atoms_from = 0
676 colvar%population_param%use_kinds_from = .true.
679 CALL uppercase(colvar%population_param%c_kinds_from(k))
689 CALL reallocate(colvar%population_param%i_at_to, 1, ndim +
SIZE(iatms))
690 colvar%population_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
691 ndim = ndim +
SIZE(iatms)
693 colvar%population_param%n_atoms_to = ndim
694 colvar%population_param%use_kinds_to = .false.
701 CALL reallocate(colvar%population_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
702 colvar%population_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
703 ndim = ndim +
SIZE(c_kinds)
705 colvar%population_param%n_atoms_to = 0
706 colvar%population_param%use_kinds_to = .true.
709 CALL uppercase(colvar%population_param%c_kinds_to(k))
718 ELSE IF (my_subsection(13))
THEN
720 wrk_section => plane_plane_angle_section
722 CALL colvar_check_points(colvar, plane_plane_angle_section)
727 cpabort(
"PLANE_PLANE_ANGLE Colvar section: Two PLANE sections must be provided!")
730 i_val=colvar%plane_plane_angle_param%plane1%type_of_def)
731 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_vec)
THEN
734 colvar%plane_plane_angle_param%plane1%normal_vec = s1
738 colvar%plane_plane_angle_param%plane1%points = iatms
743 i_val=colvar%plane_plane_angle_param%plane2%type_of_def)
744 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_vec)
THEN
747 colvar%plane_plane_angle_param%plane2%normal_vec = s1
751 colvar%plane_plane_angle_param%plane2%points = iatms
753 ELSE IF (my_subsection(14))
THEN
755 wrk_section => gyration_section
757 CALL colvar_check_points(colvar, gyration_section)
759 NULLIFY (colvar%gyration_param%i_at, colvar%gyration_param%c_kinds)
768 CALL reallocate(colvar%gyration_param%i_at, 1, ndim +
SIZE(iatms))
769 colvar%gyration_param%i_at(ndim + 1:ndim +
SIZE(iatms)) = iatms
770 ndim = ndim +
SIZE(iatms)
772 colvar%gyration_param%n_atoms = ndim
773 colvar%gyration_param%use_kinds = .false.
780 CALL reallocate(colvar%gyration_param%c_kinds, 1, ndim +
SIZE(c_kinds))
781 colvar%gyration_param%c_kinds(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
782 ndim = ndim +
SIZE(c_kinds)
784 colvar%gyration_param%n_atoms = 0
785 colvar%gyration_param%use_kinds = .true.
788 CALL uppercase(colvar%gyration_param%c_kinds(k))
791 ELSE IF (my_subsection(15))
THEN
793 wrk_section => rmsd_section
796 NULLIFY (colvar%rmsd_param%i_rmsd, colvar%rmsd_param%r_ref, colvar%rmsd_param%weights)
801 colvar%rmsd_param%nr_frames = nr_frame
803 cpassert(nr_frame >= 1 .AND. nr_frame <= 2)
804 CALL read_frames(frame_section, para_env, nr_frame, colvar%rmsd_param%r_ref, &
805 colvar%rmsd_param%n_atoms)
807 ALLOCATE (colvar%rmsd_param%weights(colvar%rmsd_param%n_atoms))
808 colvar%rmsd_param%weights = 0.0_dp
811 IF (colvar%rmsd_param%subset ==
rmsd_all)
THEN
812 ALLOCATE (colvar%rmsd_param%i_rmsd(colvar%rmsd_param%n_atoms))
813 DO i = 1, colvar%rmsd_param%n_atoms
814 colvar%rmsd_param%i_rmsd(i) = i
816 ELSE IF (colvar%rmsd_param%subset ==
rmsd_list)
THEN
824 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
825 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
826 ndim = ndim +
SIZE(iatms)
828 colvar%rmsd_param%n_atoms = ndim
830 cpabort(
"CV RMSD: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
839 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
840 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
841 ndim = ndim +
SIZE(iatms)
843 colvar%rmsd_param%n_atoms = ndim
845 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of atoms needs to be provided ")
853 CALL reallocate(weights, 1, ndim +
SIZE(wei))
854 weights(ndim + 1:ndim +
SIZE(wei)) = wei
855 ndim = ndim +
SIZE(wei)
857 IF (ndim /= colvar%rmsd_param%n_atoms) &
858 CALL cp_abort(__location__,
"CV RMSD: list of atoms and list of "// &
859 "weights need to contain same number of entries. ")
861 ii = colvar%rmsd_param%i_rmsd(i)
862 colvar%rmsd_param%weights(ii) = weights(i)
866 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of weights need to be provided. ")
870 cpabort(
"CV RMSD: unknown SUBSET_TYPE.")
874 ELSE IF (my_subsection(17))
THEN
876 wrk_section => xyz_diag_section
878 CALL colvar_check_points(colvar, wrk_section)
882 CALL section_vals_val_get(wrk_section,
"ABSOLUTE_POSITION", l_val=colvar%xyz_diag_param%use_absolute_position)
883 colvar%xyz_diag_param%i_atom = iatm
884 colvar%xyz_diag_param%component = icomponent
885 ELSE IF (my_subsection(18))
THEN
887 wrk_section => xyz_outerdiag_section
889 CALL colvar_check_points(colvar, wrk_section)
891 colvar%xyz_outerdiag_param%i_atoms = iatms
893 colvar%xyz_outerdiag_param%components(1) = icomponent
895 colvar%xyz_outerdiag_param%components(2) = icomponent
897 ELSE IF (my_subsection(19))
THEN
899 wrk_section => u_section
902 CALL section_vals_get(colvar%u_param%mixed_energy_section, explicit=use_mixed_energy)
903 IF (.NOT. use_mixed_energy)
NULLIFY (colvar%u_param%mixed_energy_section)
904 ELSE IF (my_subsection(20))
THEN
906 wrk_section => wc_section
908 CALL colvar_check_points(colvar, wc_section)
912 colvar%Wc%ids = iatms
913 ELSE IF (my_subsection(21))
THEN
915 wrk_section => hbp_section
917 CALL colvar_check_points(colvar, hbp_section)
923 ALLOCATE (colvar%HBP%ids(colvar%HBP%nPoints, 3))
924 ALLOCATE (colvar%HBP%ewc(colvar%HBP%nPoints))
925 DO i = 1, colvar%HBP%nPoints
927 colvar%HBP%ids(i, :) = iatms
929 ELSE IF (my_subsection(22))
THEN
933 colvar%ring_puckering_param%nring =
SIZE(iatms)
934 ALLOCATE (colvar%ring_puckering_param%atoms(
SIZE(iatms)))
935 colvar%ring_puckering_param%atoms = iatms
937 i_val=colvar%ring_puckering_param%iq)
939 ndim = colvar%ring_puckering_param%nring
941 cpabort(
"CV Ring Puckering: Ring size has to be 4 or larger. ")
942 ii = colvar%ring_puckering_param%iq
943 IF (abs(ii) == 1 .OR. ii < -(ndim - 1)/2 .OR. ii > ndim/2) &
944 cpabort(
"CV Ring Puckering: Invalid coordinate number.")
945 ELSE IF (my_subsection(23))
THEN
947 wrk_section => mindist_section
949 CALL colvar_check_points(colvar, mindist_section)
950 NULLIFY (colvar%mindist_param%i_dist_from, colvar%mindist_param%i_coord_from, &
951 colvar%mindist_param%k_coord_from, colvar%mindist_param%i_coord_to, &
952 colvar%mindist_param%k_coord_to)
954 colvar%mindist_param%n_dist_from =
SIZE(iatms)
955 ALLOCATE (colvar%mindist_param%i_dist_from(
SIZE(iatms)))
956 colvar%mindist_param%i_dist_from = iatms
963 CALL reallocate(colvar%mindist_param%i_coord_from, 1, ndim +
SIZE(iatms))
964 colvar%mindist_param%i_coord_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
965 ndim = ndim +
SIZE(iatms)
967 colvar%mindist_param%n_coord_from = ndim
968 colvar%mindist_param%use_kinds_from = .false.
975 CALL reallocate(colvar%mindist_param%k_coord_from, 1, ndim +
SIZE(c_kinds))
976 colvar%mindist_param%k_coord_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
977 ndim = ndim +
SIZE(c_kinds)
979 colvar%mindist_param%n_coord_from = 0
980 colvar%mindist_param%use_kinds_from = .true.
983 CALL uppercase(colvar%mindist_param%k_coord_from(k))
993 CALL reallocate(colvar%mindist_param%i_coord_to, 1, ndim +
SIZE(iatms))
994 colvar%mindist_param%i_coord_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
995 ndim = ndim +
SIZE(iatms)
997 colvar%mindist_param%n_coord_to = ndim
998 colvar%mindist_param%use_kinds_to = .false.
1005 CALL reallocate(colvar%mindist_param%k_coord_to, 1, ndim +
SIZE(c_kinds))
1006 colvar%mindist_param%k_coord_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
1007 ndim = ndim +
SIZE(c_kinds)
1009 colvar%mindist_param%n_coord_to = 0
1010 colvar%mindist_param%use_kinds_to = .true.
1013 CALL uppercase(colvar%mindist_param%k_coord_to(k))
1022 ELSE IF (my_subsection(24))
THEN
1025 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_water)
1026 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_acid)
1027 NULLIFY (colvar%acid_hyd_dist_param%i_hydrogens)
1029 colvar%acid_hyd_dist_param%n_oxygens_water, &
1030 colvar%acid_hyd_dist_param%n_oxygens_acid, &
1031 colvar%acid_hyd_dist_param%n_hydrogens, &
1032 colvar%acid_hyd_dist_param%i_oxygens_water, &
1033 colvar%acid_hyd_dist_param%i_oxygens_acid, &
1034 colvar%acid_hyd_dist_param%i_hydrogens)
1035 ELSE IF (my_subsection(25))
THEN
1038 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_water)
1039 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_acid)
1040 NULLIFY (colvar%acid_hyd_shell_param%i_hydrogens)
1042 colvar%acid_hyd_shell_param%n_oxygens_water, &
1043 colvar%acid_hyd_shell_param%n_oxygens_acid, &
1044 colvar%acid_hyd_shell_param%n_hydrogens, &
1045 colvar%acid_hyd_shell_param%i_oxygens_water, &
1046 colvar%acid_hyd_shell_param%i_oxygens_acid, &
1047 colvar%acid_hyd_shell_param%i_hydrogens)
1048 ELSE IF (my_subsection(26))
THEN
1051 NULLIFY (colvar%hydronium_dist_param%i_oxygens)
1052 NULLIFY (colvar%hydronium_dist_param%i_hydrogens)
1054 colvar%hydronium_dist_param%n_oxygens, &
1055 colvar%hydronium_dist_param%n_hydrogens, &
1056 colvar%hydronium_dist_param%i_oxygens, &
1057 colvar%hydronium_dist_param%i_hydrogens)
1062 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
1065 IF (colvar%use_points) tag =
"POINTS:"
1068 WRITE (iw,
'( A )')
' '// &
1069 '----------------------------------------------------------------------'
1070 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
1073 SELECT CASE (colvar%type_id)
1075 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| ANGLE >>> '//tag, &
1076 colvar%angle_param%i_at_angle
1078 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| DISTANCE DIFFERENCE >>> '//tag, &
1079 colvar%dfunct_param%i_at_dfunct
1081 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE DISTANCE - PLANE >>> '//tag, &
1082 colvar%plane_distance_param%plane
1083 WRITE (iw,
'( A,T73,1I8)')
' COLVARS| PLANE DISTANCE - POINT >>> '//tag, &
1084 colvar%plane_distance_param%point
1086 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
1087 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1088 colvar%plane_plane_angle_param%plane1%points
1090 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1091 colvar%plane_plane_angle_param%plane1%normal_vec
1094 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
1095 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1096 colvar%plane_plane_angle_param%plane2%points
1098 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1099 colvar%plane_plane_angle_param%plane2%normal_vec
1102 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| TORSION >>> '//tag, &
1103 colvar%torsion_param%i_at_tors
1105 WRITE (iw,
'( A,T65,2I8)')
' COLVARS| BOND >>> '//tag, &
1106 colvar%dist_param%i_at, colvar%dist_param%j_at
1108 IF (colvar%coord_param%do_chain)
THEN
1109 WRITE (iw,
'( A)')
' COLVARS| COORDINATION CHAIN FC(from->to)*FC(to->to_B)>> '
1111 IF (colvar%coord_param%use_kinds_from)
THEN
1112 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> FROM KINDS', &
1113 adjustr(colvar%coord_param%c_kinds_from(kk) (1:10)), &
1114 kk=1,
SIZE(colvar%coord_param%c_kinds_from))
1116 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> FROM '//tag, &
1117 colvar%coord_param%i_at_from(kk), &
1118 kk=1,
SIZE(colvar%coord_param%i_at_from))
1120 IF (colvar%coord_param%use_kinds_to)
THEN
1121 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS', &
1122 adjustr(colvar%coord_param%c_kinds_to(kk) (1:10)), &
1123 kk=1,
SIZE(colvar%coord_param%c_kinds_to))
1125 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag, &
1126 colvar%coord_param%i_at_to(kk), &
1127 kk=1,
SIZE(colvar%coord_param%i_at_to))
1129 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%coord_param%r_0
1130 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%coord_param%nncrd
1131 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%coord_param%ndcrd
1132 IF (colvar%coord_param%do_chain)
THEN
1133 IF (colvar%coord_param%use_kinds_to_b)
THEN
1134 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS B', &
1135 adjustr(colvar%coord_param%c_kinds_to_b(kk) (1:10)), &
1136 kk=1,
SIZE(colvar%coord_param%c_kinds_to_b))
1138 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag//
' B', &
1139 colvar%coord_param%i_at_to_b(kk), &
1140 kk=1,
SIZE(colvar%coord_param%i_at_to_b))
1142 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0 B', colvar%coord_param%r_0_b
1143 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN B', colvar%coord_param%nncrd_b
1144 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND B', colvar%coord_param%ndcrd_b
1147 IF (colvar%population_param%use_kinds_from)
THEN
1148 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> FROM KINDS', &
1149 adjustr(colvar%population_param%c_kinds_from(kk) (1:10)), &
1150 kk=1,
SIZE(colvar%population_param%c_kinds_from))
1152 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> FROM '//tag, &
1153 colvar%population_param%i_at_from(kk), &
1154 kk=1,
SIZE(colvar%population_param%i_at_from))
1156 IF (colvar%population_param%use_kinds_to)
THEN
1157 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> TO KINDS', &
1158 adjustr(colvar%population_param%c_kinds_to(kk) (1:10)), &
1159 kk=1,
SIZE(colvar%population_param%c_kinds_to))
1161 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> TO '//tag, &
1162 colvar%population_param%i_at_to(kk), &
1163 kk=1,
SIZE(colvar%population_param%i_at_to))
1165 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%population_param%r_0
1166 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%population_param%nncrd
1167 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%population_param%ndcrd
1168 WRITE (iw,
'( A,T71,I10)')
' COLVARS| N0', colvar%population_param%n0
1169 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| SIGMA', colvar%population_param%sigma
1171 IF (colvar%gyration_param%use_kinds)
THEN
1172 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| Gyration Radius >>> KINDS', &
1173 adjustr(colvar%gyration_param%c_kinds(kk) (1:10)), &
1174 kk=1,
SIZE(colvar%gyration_param%c_kinds))
1176 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Gyration Radius >>> ATOMS '//tag, &
1177 colvar%gyration_param%i_at(kk), &
1178 kk=1,
SIZE(colvar%gyration_param%i_at))
1181 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 1 >>> '//tag, &
1182 colvar%rotation_param%i_at1_bond1
1183 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 1 >>> '//tag, &
1184 colvar%rotation_param%i_at2_bond1
1185 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 2 >>> '//tag, &
1186 colvar%rotation_param%i_at1_bond2
1187 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 2 >>> '//tag, &
1188 colvar%rotation_param%i_at2_bond2
1190 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> FROM '//tag, &
1191 colvar%qparm_param%i_at_from(kk), &
1192 kk=1,
SIZE(colvar%qparm_param%i_at_from))
1193 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> TO '//tag, &
1194 colvar%qparm_param%i_at_to(kk), &
1195 kk=1,
SIZE(colvar%qparm_param%i_at_to))
1196 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RCUT', colvar%qparm_param%rcut
1197 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RSTART', colvar%qparm_param%rstart
1198 WRITE (iw,
'( A,T71,L10)')
' COLVARS| INCLUDE IMAGES', colvar%qparm_param%include_images
1200 WRITE (iw,
'( A,T71,I10)')
' COLVARS| L', colvar%qparm_param%l
1202 WRITE (iw,
'( A)')
' COLVARS| COMBINING FUNCTION : '// &
1203 trim(colvar%combine_cvs_param%function)
1204 WRITE (iw,
'( A)', advance=
"NO")
' COLVARS| VARIABLES : '
1205 DO i = 1,
SIZE(colvar%combine_cvs_param%variables)
1206 WRITE (iw,
'( A)', advance=
"NO") &
1207 trim(colvar%combine_cvs_param%variables(i))//
" "
1210 WRITE (iw,
'( A)')
' COLVARS| DEFINED PARAMETERS [label] [value]:'
1211 DO i = 1,
SIZE(colvar%combine_cvs_param%c_parameters)
1212 WRITE (iw,
'( A,A7,F9.3)')
' ', &
1213 trim(colvar%combine_cvs_param%c_parameters(i)), colvar%combine_cvs_param%v_parameters(i)
1215 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| ERROR ON DERIVATIVE EVALUATION', &
1216 colvar%combine_cvs_param%lerr
1217 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| DX', &
1218 colvar%combine_cvs_param%dx
1220 cpwarn(
"Description header for REACTION_PATH COLVAR missing!!")
1222 cpwarn(
"Description header for REACTION_PATH COLVAR missing!!")
1224 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_shell_param%poh
1225 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_shell_param%qoh
1226 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%hydronium_shell_param%poo
1227 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%hydronium_shell_param%qoo
1228 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%hydronium_shell_param%roo
1229 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_shell_param%roh
1230 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_shell_param%nh
1231 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%hydronium_shell_param%lambda
1233 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_dist_param%poh
1234 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_dist_param%qoh
1235 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_dist_param%roh
1236 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%hydronium_dist_param%pm
1237 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%hydronium_dist_param%qm
1238 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_dist_param%nh
1239 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PF', colvar%hydronium_dist_param%pf
1240 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QF', colvar%hydronium_dist_param%qf
1241 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NN', colvar%hydronium_dist_param%nn
1243 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_dist_param%paoh
1244 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_dist_param%qaoh
1245 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_dist_param%pwoh
1246 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_dist_param%qwoh
1247 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_dist_param%pcut
1248 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_dist_param%qcut
1249 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_dist_param%raoh
1250 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_dist_param%rwoh
1251 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_dist_param%nc
1252 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_dist_param%lambda
1254 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_shell_param%paoh
1255 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_shell_param%qaoh
1256 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_shell_param%pwoh
1257 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_shell_param%qwoh
1258 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%acid_hyd_shell_param%poo
1259 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%acid_hyd_shell_param%qoo
1260 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%acid_hyd_shell_param%pm
1261 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%acid_hyd_shell_param%qm
1262 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_shell_param%pcut
1263 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_shell_param%qcut
1264 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_shell_param%raoh
1265 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_shell_param%rwoh
1266 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%acid_hyd_shell_param%roo
1267 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%acid_hyd_shell_param%nh
1268 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_shell_param%nc
1269 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_shell_param%lambda
1271 cpwarn(
"Description header for RMSD COLVAR missing!!")
1273 NULLIFY (section, keyword, enum)
1277 tag_comp =
enum_i2c(enum, colvar%xyz_diag_param%component)
1280 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| POSITION ('//trim(tag_comp) &
1281 //
') >>> '//tag, colvar%xyz_diag_param%i_atom
1283 NULLIFY (section, keyword, enum)
1287 tag_comp1 =
enum_i2c(enum, colvar%xyz_outerdiag_param%components(1))
1290 tag_comp2 =
enum_i2c(enum, colvar%xyz_outerdiag_param%components(2))
1293 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| CROSS TERM POSITION ('//trim(tag_comp1) &
1294 //
" * "//trim(tag_comp2)//
') >>> '//tag, colvar%xyz_outerdiag_param%i_atoms
1296 WRITE (iw,
'( A,T77,A4)')
' COLVARS| ENERGY >>> '//tag,
'all!'
1298 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| Wc >>> RCUT: ', &
1300 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| Wc >>> '//tag, &
1303 WRITE (iw,
'( A,T57,I8)')
' COLVARS| HBP >>> NPOINTS', &
1305 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1307 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1309 DO i = 1, colvar%HBP%nPoints
1310 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| HBP >>> '//tag, &
1311 colvar%HBP%ids(i, :)
1314 WRITE (iw,
'( A,T57,I8)')
' COLVARS| Ring Puckering >>> ring size', &
1315 colvar%ring_puckering_param%nring
1316 IF (colvar%ring_puckering_param%iq == 0)
THEN
1317 WRITE (iw,
'( A,T40,A)')
' COLVARS| Ring Puckering >>> coordinate', &
1318 ' Total Puckering Amplitude'
1319 ELSEIF (colvar%ring_puckering_param%iq > 0)
THEN
1320 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1321 ' Puckering Amplitude', &
1322 colvar%ring_puckering_param%iq
1324 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1325 ' Puckering Angle', &
1326 colvar%ring_puckering_param%iq
1329 WRITE (iw,
'( A)')
' COLVARS| CONDITIONED DISTANCE>> '
1330 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DISTANCE >>> DISTANCE FROM '//tag, &
1331 colvar%mindist_param%i_dist_from(kk), &
1332 kk=1,
SIZE(colvar%mindist_param%i_dist_from))
1333 IF (colvar%mindist_param%use_kinds_from)
THEN
1334 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM KINDS ', &
1335 adjustr(colvar%mindist_param%k_coord_from(kk) (1:10)), &
1336 kk=1,
SIZE(colvar%mindist_param%k_coord_from))
1338 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM '//tag, &
1339 colvar%mindist_param%i_coord_from(kk), &
1340 kk=1,
SIZE(colvar%mindist_param%i_coord_from))
1342 IF (colvar%mindist_param%use_kinds_to)
THEN
1343 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO KINDS ', &
1344 adjustr(colvar%mindist_param%k_coord_to(kk) (1:10)), &
1345 kk=1,
SIZE(colvar%mindist_param%k_coord_to))
1347 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO '//tag, &
1348 colvar%mindist_param%i_coord_to(kk), &
1349 kk=1,
SIZE(colvar%mindist_param%i_coord_to))
1351 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%mindist_param%r_cut
1352 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%mindist_param%p_exp
1353 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%mindist_param%q_exp
1354 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%mindist_param%lambda
1357 IF (colvar%use_points)
THEN
1358 WRITE (iw,
'( A)')
' COLVARS| INFORMATION ON DEFINED GEOMETRICAL POINTS'
1359 DO kk = 1,
SIZE(colvar%points)
1362 tmpstr2 = cp_to_string(kk)
1363 WRITE (iw,
'( A)')
' COLVARS| POINT Nr.'//trim(tmpstr2)//
' OF TYPE: '//trim(tmpstr)
1364 IF (
ASSOCIATED(colvar%points(kk)%atoms))
THEN
1365 WRITE (iw,
'( A)')
' COLVARS| ATOMS BUILDING THE GEOMETRICAL POINT'
1366 WRITE (iw,
'( A, I10)') (
' COLVARS| ATOM:', colvar%points(kk)%atoms(k), k=1,
SIZE(colvar%points(kk)%atoms))
1368 WRITE (iw,
'( A,4X,3F12.6)')
' COLVARS| XYZ POSITION OF FIXED POINT:', colvar%points(kk)%r
1374 WRITE (iw,
'( A )')
' '// &
1375 '----------------------------------------------------------------------'
1377 WRITE (iw,
'( A )')
' '// &
1378 '**********************************************************************'
1382 "PRINT%PROGRAM_RUN_INFO")
1383 CALL timestop(handle)
1397 SUBROUTINE read_hydronium_colvars(section, colvar, colvar_id, n_oxygens, n_hydrogens, &
1398 i_oxygens, i_hydrogens)
1399 TYPE(section_vals_type),
POINTER :: section
1400 TYPE(colvar_type),
POINTER :: colvar
1401 INTEGER,
INTENT(IN) :: colvar_id
1402 INTEGER,
INTENT(OUT) :: n_oxygens, n_hydrogens
1403 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens, i_hydrogens
1405 INTEGER :: k, n_var, ndim
1406 INTEGER,
DIMENSION(:),
POINTER :: iatms
1414 CALL reallocate(i_oxygens, 1, ndim +
SIZE(iatms))
1415 i_oxygens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1416 ndim = ndim +
SIZE(iatms)
1424 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1425 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1426 ndim = ndim +
SIZE(iatms)
1430 SELECT CASE (colvar_id)
1455 END SUBROUTINE read_hydronium_colvars
1471 SUBROUTINE read_acid_hydronium_colvars(section, colvar, colvar_id, n_oxygens_water, &
1472 n_oxygens_acid, n_hydrogens, i_oxygens_water, &
1473 i_oxygens_acid, i_hydrogens)
1474 TYPE(section_vals_type),
POINTER :: section
1475 TYPE(colvar_type),
POINTER :: colvar
1476 INTEGER,
INTENT(IN) :: colvar_id
1477 INTEGER,
INTENT(OUT) :: n_oxygens_water, n_oxygens_acid, &
1479 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens_water, i_oxygens_acid, &
1482 INTEGER :: k, n_var, ndim
1483 INTEGER,
DIMENSION(:),
POINTER :: iatms
1491 CALL reallocate(i_oxygens_water, 1, ndim +
SIZE(iatms))
1492 i_oxygens_water(ndim + 1:ndim +
SIZE(iatms)) = iatms
1493 ndim = ndim +
SIZE(iatms)
1495 n_oxygens_water = ndim
1501 CALL reallocate(i_oxygens_acid, 1, ndim +
SIZE(iatms))
1502 i_oxygens_acid(ndim + 1:ndim +
SIZE(iatms)) = iatms
1503 ndim = ndim +
SIZE(iatms)
1505 n_oxygens_acid = ndim
1511 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1512 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1513 ndim = ndim +
SIZE(iatms)
1517 SELECT CASE (colvar_id)
1548 END SUBROUTINE read_acid_hydronium_colvars
1556 SUBROUTINE colvar_check_points(colvar, section)
1557 TYPE(colvar_type),
POINTER :: colvar
1558 TYPE(section_vals_type),
POINTER :: section
1560 INTEGER :: i, irep, natoms, npoints, nrep, nweights
1561 INTEGER,
DIMENSION(:),
POINTER :: atoms
1563 REAL(kind=
dp),
DIMENSION(:),
POINTER :: r, weights
1564 TYPE(section_vals_type),
POINTER :: point_sections
1566 NULLIFY (point_sections)
1569 cpassert(
ASSOCIATED(colvar))
1573 colvar%use_points = .true.
1575 ALLOCATE (colvar%points(npoints))
1580 NULLIFY (colvar%points(i)%atoms)
1581 NULLIFY (colvar%points(i)%weights)
1582 CALL section_vals_val_get(point_sections,
"TYPE", i_rep_section=i, i_val=colvar%points(i)%type_id)
1583 SELECT CASE (colvar%points(i)%type_id)
1586 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, n_rep_val=nrep, i_vals=atoms)
1588 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1589 natoms = natoms +
SIZE(atoms)
1591 ALLOCATE (colvar%points(i)%atoms(natoms))
1594 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1595 colvar%points(i)%atoms(natoms + 1:) = atoms(:)
1596 natoms = natoms +
SIZE(atoms)
1599 ALLOCATE (colvar%points(i)%weights(natoms))
1600 colvar%points(i)%weights = 1.0_dp/real(natoms, kind=
dp)
1606 colvar%points(i)%weights(nweights + 1:) = weights(:)
1607 nweights = nweights +
SIZE(weights)
1609 cpassert(natoms == nweights)
1614 colvar%points(i)%r = r
1618 END SUBROUTINE colvar_check_points
1631 TYPE(colvar_type),
POINTER :: colvar
1632 TYPE(cell_type),
POINTER :: cell
1633 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
1634 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
1636 TYPE(fixd_constraint_type),
DIMENSION(:), &
1637 OPTIONAL,
POINTER :: fixd_list
1640 LOGICAL :: colvar_ok
1642 colvar_ok =
ASSOCIATED(colvar)
1645 IF (
PRESENT(pos))
THEN
1646 DO i = 1,
SIZE(colvar%i_atom)
1647 j = colvar%i_atom(i)
1648 particles(j)%r = pos(:, j)
1652 colvar%dsdr = 0.0_dp
1653 SELECT CASE (colvar%type_id)
1655 CALL dist_colvar(colvar, cell, particles=particles)
1657 CALL coord_colvar(colvar, cell, particles=particles)
1659 CALL population_colvar(colvar, cell, particles=particles)
1661 CALL gyration_radius_colvar(colvar, cell, particles=particles)
1663 CALL torsion_colvar(colvar, cell, particles=particles)
1665 CALL angle_colvar(colvar, cell, particles=particles)
1667 CALL dfunct_colvar(colvar, cell, particles=particles)
1669 CALL plane_distance_colvar(colvar, cell, particles=particles)
1671 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1673 CALL rotation_colvar(colvar, cell, particles=particles)
1675 CALL qparm_colvar(colvar, cell, particles=particles)
1677 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1679 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1681 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1683 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1685 CALL rmsd_colvar(colvar, particles=particles)
1687 CALL reaction_path_colvar(colvar, cell, particles=particles)
1689 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1691 CALL combine_colvar(colvar, cell, particles=particles)
1693 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1695 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1697 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1699 CALL mindist_colvar(colvar, cell, particles=particles)
1701 cpabort(
"need force_env!")
1704 CALL wc_colvar(colvar, cell, particles=particles)
1707 CALL hbp_colvar(colvar, cell, particles=particles)
1727 TYPE(force_env_type),
POINTER :: force_env
1729 LOGICAL :: colvar_ok
1730 TYPE(cell_type),
POINTER :: cell
1731 TYPE(colvar_type),
POINTER :: colvar
1732 TYPE(cp_subsys_type),
POINTER :: subsys
1733 TYPE(qs_environment_type),
POINTER :: qs_env
1735 NULLIFY (subsys, cell, colvar, qs_env)
1736 CALL force_env_get(force_env, subsys=subsys, cell=cell, qs_env=qs_env)
1737 colvar_ok =
ASSOCIATED(subsys%colvar_p)
1740 colvar => subsys%colvar_p(icolvar)%colvar
1742 colvar%dsdr = 0.0_dp
1743 SELECT CASE (colvar%type_id)
1745 CALL dist_colvar(colvar, cell, subsys=subsys)
1747 CALL coord_colvar(colvar, cell, subsys=subsys)
1749 CALL population_colvar(colvar, cell, subsys=subsys)
1751 CALL gyration_radius_colvar(colvar, cell, subsys=subsys)
1753 CALL torsion_colvar(colvar, cell, subsys=subsys, no_riemann_sheet_op=.true.)
1755 CALL angle_colvar(colvar, cell, subsys=subsys)
1757 CALL dfunct_colvar(colvar, cell, subsys=subsys)
1759 CALL plane_distance_colvar(colvar, cell, subsys=subsys)
1761 CALL plane_plane_angle_colvar(colvar, cell, subsys=subsys)
1763 CALL rotation_colvar(colvar, cell, subsys=subsys)
1765 CALL qparm_colvar(colvar, cell, subsys=subsys)
1767 CALL hydronium_shell_colvar(colvar, cell, subsys=subsys)
1769 CALL hydronium_dist_colvar(colvar, cell, subsys=subsys)
1771 CALL acid_hyd_dist_colvar(colvar, cell, subsys=subsys)
1773 CALL acid_hyd_shell_colvar(colvar, cell, subsys=subsys)
1775 CALL rmsd_colvar(colvar, subsys=subsys)
1777 CALL reaction_path_colvar(colvar, cell, subsys=subsys)
1779 CALL distance_from_path_colvar(colvar, cell, subsys=subsys)
1781 CALL combine_colvar(colvar, cell, subsys=subsys)
1783 CALL xyz_diag_colvar(colvar, cell, subsys=subsys)
1785 CALL xyz_outerdiag_colvar(colvar, cell, subsys=subsys)
1787 CALL u_colvar(colvar, force_env=force_env)
1789 CALL wc_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1791 CALL hbp_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1793 CALL ring_puckering_colvar(colvar, cell, subsys=subsys)
1795 CALL mindist_colvar(colvar, cell, subsys=subsys)
1811 SUBROUTINE colvar_recursive_eval(colvar, cell, particles)
1812 TYPE(colvar_type),
POINTER :: colvar
1813 TYPE(cell_type),
POINTER :: cell
1814 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
1818 colvar%dsdr = 0.0_dp
1819 SELECT CASE (colvar%type_id)
1821 CALL dist_colvar(colvar, cell, particles=particles)
1823 CALL coord_colvar(colvar, cell, particles=particles)
1825 CALL torsion_colvar(colvar, cell, particles=particles)
1827 CALL angle_colvar(colvar, cell, particles=particles)
1829 CALL dfunct_colvar(colvar, cell, particles=particles)
1831 CALL plane_distance_colvar(colvar, cell, particles=particles)
1833 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1835 CALL rotation_colvar(colvar, cell, particles=particles)
1837 CALL qparm_colvar(colvar, cell, particles=particles)
1839 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1841 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1843 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1845 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1847 CALL rmsd_colvar(colvar, particles=particles)
1849 CALL reaction_path_colvar(colvar, cell, particles=particles)
1851 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1853 CALL combine_colvar(colvar, cell, particles=particles)
1855 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1857 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1859 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1861 CALL mindist_colvar(colvar, cell, particles=particles)
1863 cpabort(
"need force_env!")
1865 CALL wc_colvar(colvar, cell, particles=particles)
1867 CALL hbp_colvar(colvar, cell, particles=particles)
1871 END SUBROUTINE colvar_recursive_eval
1881 SUBROUTINE get_coordinates(colvar, i, ri, my_particles)
1882 TYPE(colvar_type),
POINTER :: colvar
1883 INTEGER,
INTENT(IN) :: i
1884 REAL(kind=
dp),
DIMENSION(3),
INTENT(OUT) :: ri
1885 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
1887 IF (colvar%use_points)
THEN
1890 ri(:) = my_particles(i)%r(:)
1893 END SUBROUTINE get_coordinates
1903 SUBROUTINE get_mass(colvar, i, mi, my_particles)
1904 TYPE(colvar_type),
POINTER :: colvar
1905 INTEGER,
INTENT(IN) :: i
1906 REAL(kind=
dp),
INTENT(OUT) :: mi
1907 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
1909 IF (colvar%use_points)
THEN
1912 mi = my_particles(i)%atomic_kind%mass
1915 END SUBROUTINE get_mass
1924 SUBROUTINE put_derivative(colvar, i, fi)
1925 TYPE(colvar_type),
POINTER :: colvar
1926 INTEGER,
INTENT(IN) :: i
1927 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: fi
1929 IF (colvar%use_points)
THEN
1932 colvar%dsdr(:, i) = colvar%dsdr(:, i) + fi
1935 END SUBROUTINE put_derivative
1945 SUBROUTINE xyz_diag_colvar(colvar, cell, subsys, particles)
1946 TYPE(colvar_type),
POINTER :: colvar
1947 TYPE(cell_type),
POINTER :: cell
1948 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
1949 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
1950 POINTER :: particles
1953 REAL(
dp) :: fi(3), r, r0(3), ss(3), xi(3), xpi(3)
1954 TYPE(particle_list_type),
POINTER :: particles_i
1955 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
1957 NULLIFY (particles_i)
1960 IF (
PRESENT(particles))
THEN
1961 my_particles => particles
1963 cpassert(
PRESENT(subsys))
1965 my_particles => particles_i%els
1967 i = colvar%xyz_diag_param%i_atom
1969 CALL get_coordinates(colvar, i, xpi, my_particles)
1972 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
1973 IF (all(colvar%xyz_diag_param%r0 == huge(0.0_dp)))
THEN
1974 colvar%xyz_diag_param%r0 = xpi
1976 r0 = colvar%xyz_diag_param%r0
1981 IF (colvar%xyz_diag_param%use_pbc)
THEN
1982 ss = matmul(cell%h_inv, xpi - r0)
1984 xi = matmul(cell%hmat, ss)
1989 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
1990 SELECT CASE (colvar%xyz_diag_param%component)
2010 r = xi(1)**2 + xi(2)**2 + xi(3)**2
2013 SELECT CASE (colvar%xyz_diag_param%component)
2037 CALL put_derivative(colvar, 1, fi)
2039 END SUBROUTINE xyz_diag_colvar
2049 SUBROUTINE xyz_outerdiag_colvar(colvar, cell, subsys, particles)
2050 TYPE(colvar_type),
POINTER :: colvar
2051 TYPE(cell_type),
POINTER :: cell
2052 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2053 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2054 POINTER :: particles
2057 REAL(
dp) :: fi(3, 2), r, r0(3), ss(3), xi(3, 2), &
2059 TYPE(particle_list_type),
POINTER :: particles_i
2060 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2062 NULLIFY (particles_i)
2065 IF (
PRESENT(particles))
THEN
2066 my_particles => particles
2068 cpassert(
PRESENT(subsys))
2070 my_particles => particles_i%els
2073 i = colvar%xyz_outerdiag_param%i_atoms(k)
2075 CALL get_coordinates(colvar, i, xpi, my_particles)
2076 r0 = colvar%xyz_outerdiag_param%r0(:, k)
2077 IF (all(colvar%xyz_outerdiag_param%r0(:, k) == huge(0.0_dp))) r0 = xpi
2079 IF (colvar%xyz_outerdiag_param%use_pbc)
THEN
2080 ss = matmul(cell%h_inv, xpi - r0)
2082 xi(:, k) = matmul(cell%hmat, ss)
2087 SELECT CASE (colvar%xyz_outerdiag_param%components(k))
2112 IF (xi(l, 1) /= 0.0_dp) fi(l, 1) = fi(l, 1) + xi(i, 2)
2113 r = r + xi(l, 1)*xi(i, 2)
2115 IF (xi(i, 2) /= 0.0_dp) fi(i, 2) = sum(xi(:, 1))
2119 CALL put_derivative(colvar, 1, fi(:, 1))
2120 CALL put_derivative(colvar, 2, fi(:, 2))
2122 END SUBROUTINE xyz_outerdiag_colvar
2132 SUBROUTINE u_colvar(colvar, force_env)
2133 TYPE(colvar_type),
POINTER :: colvar
2134 TYPE(force_env_type),
OPTIONAL,
POINTER :: force_env
2136 CHARACTER(LEN=default_path_length) :: coupling_function
2137 CHARACTER(LEN=default_string_length) :: def_error, this_error
2138 CHARACTER(LEN=default_string_length), &
2139 DIMENSION(:),
POINTER :: parameters
2140 INTEGER :: iatom, iforce_eval, iparticle, &
2141 jparticle, natom, natom_iforce, &
2143 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms, map_index
2144 REAL(
dp) :: dedf, dx, err, fi(3), lerr, &
2146 REAL(kind=
dp),
DIMENSION(:),
POINTER :: values
2147 TYPE(cp_subsys_p_type),
DIMENSION(:),
POINTER :: subsystems
2148 TYPE(cp_subsys_type),
POINTER :: subsys_main
2149 TYPE(mixed_force_type),
DIMENSION(:),
POINTER :: global_forces
2150 TYPE(particle_list_p_type),
DIMENSION(:),
POINTER :: particles
2151 TYPE(particle_list_type),
POINTER :: particles_main
2152 TYPE(section_vals_type),
POINTER :: force_env_section, mapping_section, &
2155 IF (
PRESENT(force_env))
THEN
2156 NULLIFY (particles_main, subsys_main)
2158 CALL cp_subsys_get(subsys=subsys_main, particles=particles_main)
2159 natom =
SIZE(particles_main%els)
2160 colvar%n_atom_s = natom
2161 colvar%u_param%natom = natom
2162 CALL reallocate(colvar%i_atom, 1, natom)
2163 CALL reallocate(colvar%dsdr, 1, 3, 1, natom)
2165 colvar%i_atom(iatom) = iatom
2168 IF (.NOT.
ASSOCIATED(colvar%u_param%mixed_energy_section))
THEN
2169 CALL force_env_get(force_env, potential_energy=potential_energy)
2170 colvar%ss = potential_energy
2174 fi(:) = -particles_main%els(iatom)%f
2175 CALL put_derivative(colvar, iatom, fi)
2179 CALL cp_abort(__location__, &
2180 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
2181 ' A combination of mixed force_eval energies has been requested as '// &
2182 ' collective variable, but the MIXED env is not in use! Aborting.')
2183 CALL force_env_get(force_env, force_env_section=force_env_section)
2185 NULLIFY (values, parameters, subsystems, particles, global_forces, map_index, glob_natoms)
2186 nforce_eval =
SIZE(force_env%sub_force_env)
2187 ALLOCATE (glob_natoms(nforce_eval))
2188 ALLOCATE (subsystems(nforce_eval))
2189 ALLOCATE (particles(nforce_eval))
2191 ALLOCATE (global_forces(nforce_eval))
2194 DO iforce_eval = 1, nforce_eval
2195 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
2196 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
2198 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
2199 subsys=subsystems(iforce_eval)%subsys)
2202 particles=particles(iforce_eval)%list)
2205 natom_iforce =
SIZE(particles(iforce_eval)%list%els)
2208 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2209 glob_natoms(iforce_eval) = natom_iforce
2214 CALL force_env%para_env%sync()
2215 CALL force_env%para_env%sum(glob_natoms)
2218 DO iforce_eval = 1, nforce_eval
2219 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
2220 global_forces(iforce_eval)%forces = 0.0_dp
2221 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env))
THEN
2222 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2224 DO iparticle = 1, glob_natoms(iforce_eval)
2225 global_forces(iforce_eval)%forces(:, iparticle) = &
2226 particles(iforce_eval)%list%els(iparticle)%f
2230 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
2233 wrk_section => colvar%u_param%mixed_energy_section
2235 CALL get_generic_info(wrk_section,
"ENERGY_FUNCTION", coupling_function, parameters, &
2236 values, force_env%mixed_env%energies)
2238 CALL parsef(1, trim(coupling_function), parameters)
2240 colvar%ss =
evalf(1, values)
2243 DO iforce_eval = 1, nforce_eval
2246 dedf =
evalfd(1, iforce_eval, values, dx, err)
2247 IF (abs(err) > lerr)
THEN
2248 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
2249 WRITE (def_error,
"(A,G12.6,A)")
"(", lerr,
")"
2252 CALL cp_warn(__location__, &
2253 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
2254 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
2255 trim(def_error)//
' .')
2260 nforce_eval, map_index)
2263 DO iparticle = 1, glob_natoms(iforce_eval)
2264 jparticle = map_index(iparticle)
2265 fi = -dedf*global_forces(iforce_eval)%forces(:, iparticle)
2266 CALL put_derivative(colvar, jparticle, fi)
2269 IF (
ASSOCIATED(map_index))
THEN
2270 DEALLOCATE (map_index)
2274 DO iforce_eval = 1, nforce_eval
2275 DEALLOCATE (global_forces(iforce_eval)%forces)
2277 DEALLOCATE (glob_natoms)
2279 DEALLOCATE (parameters)
2280 DEALLOCATE (global_forces)
2281 DEALLOCATE (subsystems)
2282 DEALLOCATE (particles)
2285 cpabort(
"need force_env!")
2287 END SUBROUTINE u_colvar
2297 SUBROUTINE plane_distance_colvar(colvar, cell, subsys, particles)
2299 TYPE(colvar_type),
POINTER :: colvar
2300 TYPE(cell_type),
POINTER :: cell
2301 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2302 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2303 POINTER :: particles
2305 INTEGER :: i, j, k, l
2306 REAL(
dp) :: a, b, dsdxpn(3), dxpndxi(3, 3), dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), &
2307 fk(3), fl(3), r12, ri(3), rj(3), rk(3), rl(3), ss(3), xpij(3), xpkj(3), xpl(3), xpn(3)
2308 TYPE(particle_list_type),
POINTER :: particles_i
2309 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2311 NULLIFY (particles_i)
2314 IF (
PRESENT(particles))
THEN
2315 my_particles => particles
2317 cpassert(
PRESENT(subsys))
2319 my_particles => particles_i%els
2321 i = colvar%plane_distance_param%plane(1)
2322 j = colvar%plane_distance_param%plane(2)
2323 k = colvar%plane_distance_param%plane(3)
2324 l = colvar%plane_distance_param%point
2326 CALL get_coordinates(colvar, i, ri, my_particles)
2327 CALL get_coordinates(colvar, j, rj, my_particles)
2328 CALL get_coordinates(colvar, k, rk, my_particles)
2329 CALL get_coordinates(colvar, l, rl, my_particles)
2332 xpl = rl - (ri + rj + rk)/3.0_dp
2333 IF (colvar%plane_distance_param%use_pbc)
THEN
2335 ss = matmul(cell%h_inv, ri - rj)
2337 xpij = matmul(cell%hmat, ss)
2339 ss = matmul(cell%h_inv, rk - rj)
2341 xpkj = matmul(cell%hmat, ss)
2343 ss = matmul(cell%h_inv, rl - (ri + rj + rk)/3.0_dp)
2345 xpl = matmul(cell%hmat, ss)
2348 xpn(1) = xpij(2)*xpkj(3) - xpij(3)*xpkj(2)
2349 xpn(2) = xpij(3)*xpkj(1) - xpij(1)*xpkj(3)
2350 xpn(3) = xpij(1)*xpkj(2) - xpij(2)*xpkj(1)
2351 a = dot_product(xpn, xpn)
2352 b = dot_product(xpl, xpn)
2355 dsdxpn(1) = xpl(1)/r12 - b*xpn(1)/(r12*a)
2356 dsdxpn(2) = xpl(2)/r12 - b*xpn(2)/(r12*a)
2357 dsdxpn(3) = xpl(3)/r12 - b*xpn(3)/(r12*a)
2359 dxpndxi(1, 1) = 0.0_dp
2360 dxpndxi(1, 2) = 1.0_dp*xpkj(3)
2361 dxpndxi(1, 3) = -1.0_dp*xpkj(2)
2362 dxpndxi(2, 1) = -1.0_dp*xpkj(3)
2363 dxpndxi(2, 2) = 0.0_dp
2364 dxpndxi(2, 3) = 1.0_dp*xpkj(1)
2365 dxpndxi(3, 1) = 1.0_dp*xpkj(2)
2366 dxpndxi(3, 2) = -1.0_dp*xpkj(1)
2367 dxpndxi(3, 3) = 0.0_dp
2369 dxpndxj(1, 1) = 0.0_dp
2370 dxpndxj(1, 2) = -1.0_dp*xpkj(3) + xpij(3)
2371 dxpndxj(1, 3) = -1.0_dp*xpij(2) + xpkj(2)
2372 dxpndxj(2, 1) = -1.0_dp*xpij(3) + xpkj(3)
2373 dxpndxj(2, 2) = 0.0_dp
2374 dxpndxj(2, 3) = -1.0_dp*xpkj(1) + xpij(1)
2375 dxpndxj(3, 1) = -1.0_dp*xpkj(2) + xpij(2)
2376 dxpndxj(3, 2) = -1.0_dp*xpij(1) + xpkj(1)
2377 dxpndxj(3, 3) = 0.0_dp
2379 dxpndxk(1, 1) = 0.0_dp
2380 dxpndxk(1, 2) = -1.0_dp*xpij(3)
2381 dxpndxk(1, 3) = 1.0_dp*xpij(2)
2382 dxpndxk(2, 1) = 1.0_dp*xpij(3)
2383 dxpndxk(2, 2) = 0.0_dp
2384 dxpndxk(2, 3) = -1.0_dp*xpij(1)
2385 dxpndxk(3, 1) = -1.0_dp*xpij(2)
2386 dxpndxk(3, 2) = 1.0_dp*xpij(1)
2387 dxpndxk(3, 3) = 0.0_dp
2389 fi(:) = matmul(dsdxpn, dxpndxi) - xpn/(3.0_dp*r12)
2390 fj(:) = matmul(dsdxpn, dxpndxj) - xpn/(3.0_dp*r12)
2391 fk(:) = matmul(dsdxpn, dxpndxk) - xpn/(3.0_dp*r12)
2394 CALL put_derivative(colvar, 1, fi)
2395 CALL put_derivative(colvar, 2, fj)
2396 CALL put_derivative(colvar, 3, fk)
2397 CALL put_derivative(colvar, 4, fl)
2399 END SUBROUTINE plane_distance_colvar
2410 SUBROUTINE plane_plane_angle_colvar(colvar, cell, subsys, particles)
2412 TYPE(colvar_type),
POINTER :: colvar
2413 TYPE(cell_type),
POINTER :: cell
2414 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2415 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2416 POINTER :: particles
2418 INTEGER :: i1, i2, j1, j2, k1, k2, np
2420 REAL(
dp) :: a1, a2, d, dnorm_dxpn(3), dprod12_dxpn(3), dsdxpn(3), dt_dxpn(3), dxpndxi(3, 3), &
2421 dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), fk(3), fmod, norm1, norm2, prod_12, ri1(3), &
2422 ri2(3), rj1(3), rj2(3), rk1(3), rk2(3), ss(3), t, xpij1(3), xpij2(3), xpkj1(3), xpkj2(3), &
2424 TYPE(particle_list_type),
POINTER :: particles_i
2425 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2427 NULLIFY (particles_i)
2431 IF (
PRESENT(particles))
THEN
2432 my_particles => particles
2434 cpassert(
PRESENT(subsys))
2436 my_particles => particles_i%els
2440 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2441 i1 = colvar%plane_plane_angle_param%plane1%points(1)
2442 j1 = colvar%plane_plane_angle_param%plane1%points(2)
2443 k1 = colvar%plane_plane_angle_param%plane1%points(3)
2446 CALL get_coordinates(colvar, i1, ri1, my_particles)
2447 CALL get_coordinates(colvar, j1, rj1, my_particles)
2448 CALL get_coordinates(colvar, k1, rk1, my_particles)
2451 ss = matmul(cell%h_inv, ri1 - rj1)
2453 xpij1 = matmul(cell%hmat, ss)
2456 ss = matmul(cell%h_inv, rk1 - rj1)
2458 xpkj1 = matmul(cell%hmat, ss)
2461 xpn1(1) = xpij1(2)*xpkj1(3) - xpij1(3)*xpkj1(2)
2462 xpn1(2) = xpij1(3)*xpkj1(1) - xpij1(1)*xpkj1(3)
2463 xpn1(3) = xpij1(1)*xpkj1(2) - xpij1(2)*xpkj1(1)
2465 xpn1 = colvar%plane_plane_angle_param%plane1%normal_vec
2467 a1 = dot_product(xpn1, xpn1)
2469 cpassert(norm1 /= 0.0_dp)
2472 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2473 i2 = colvar%plane_plane_angle_param%plane2%points(1)
2474 j2 = colvar%plane_plane_angle_param%plane2%points(2)
2475 k2 = colvar%plane_plane_angle_param%plane2%points(3)
2478 CALL get_coordinates(colvar, i2, ri2, my_particles)
2479 CALL get_coordinates(colvar, j2, rj2, my_particles)
2480 CALL get_coordinates(colvar, k2, rk2, my_particles)
2483 ss = matmul(cell%h_inv, ri2 - rj2)
2485 xpij2 = matmul(cell%hmat, ss)
2488 ss = matmul(cell%h_inv, rk2 - rj2)
2490 xpkj2 = matmul(cell%hmat, ss)
2493 xpn2(1) = xpij2(2)*xpkj2(3) - xpij2(3)*xpkj2(2)
2494 xpn2(2) = xpij2(3)*xpkj2(1) - xpij2(1)*xpkj2(3)
2495 xpn2(3) = xpij2(1)*xpkj2(2) - xpij2(2)*xpkj2(1)
2497 xpn2 = colvar%plane_plane_angle_param%plane2%normal_vec
2499 a2 = dot_product(xpn2, xpn2)
2501 cpassert(norm2 /= 0.0_dp)
2504 prod_12 = dot_product(xpn1, xpn2)
2508 t = min(1.0_dp, abs(t))*sign(1.0_dp, t)
2511 IF ((abs(colvar%ss) .LT. tolerance_acos) .OR. (abs(colvar%ss -
pi) .LT. tolerance_acos))
THEN
2514 fmod = -1.0_dp/sin(colvar%ss)
2519 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2521 dnorm_dxpn = 1.0_dp/norm1*xpn1
2522 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm2)/d**2
2524 dsdxpn(1) = fmod*dt_dxpn(1)
2525 dsdxpn(2) = fmod*dt_dxpn(2)
2526 dsdxpn(3) = fmod*dt_dxpn(3)
2528 dxpndxi(1, 1) = 0.0_dp
2529 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2530 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2531 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2532 dxpndxi(2, 2) = 0.0_dp
2533 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2534 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2535 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2536 dxpndxi(3, 3) = 0.0_dp
2538 dxpndxj(1, 1) = 0.0_dp
2539 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2540 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2541 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2542 dxpndxj(2, 2) = 0.0_dp
2543 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2544 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2545 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2546 dxpndxj(3, 3) = 0.0_dp
2548 dxpndxk(1, 1) = 0.0_dp
2549 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2550 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2551 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2552 dxpndxk(2, 2) = 0.0_dp
2553 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2554 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2555 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2556 dxpndxk(3, 3) = 0.0_dp
2558 fi = matmul(dsdxpn, dxpndxi)
2559 fj = matmul(dsdxpn, dxpndxj)
2560 fk = matmul(dsdxpn, dxpndxk)
2563 CALL put_derivative(colvar, np + 1, fi)
2564 CALL put_derivative(colvar, np + 2, fj)
2565 CALL put_derivative(colvar, np + 3, fk)
2570 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2572 dnorm_dxpn = 1.0_dp/norm2*xpn2
2573 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm1)/d**2
2575 dsdxpn(1) = fmod*dt_dxpn(1)
2576 dsdxpn(2) = fmod*dt_dxpn(2)
2577 dsdxpn(3) = fmod*dt_dxpn(3)
2579 dxpndxi(1, 1) = 0.0_dp
2580 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2581 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2582 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2583 dxpndxi(2, 2) = 0.0_dp
2584 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2585 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2586 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2587 dxpndxi(3, 3) = 0.0_dp
2589 dxpndxj(1, 1) = 0.0_dp
2590 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2591 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2592 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2593 dxpndxj(2, 2) = 0.0_dp
2594 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2595 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2596 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2597 dxpndxj(3, 3) = 0.0_dp
2599 dxpndxk(1, 1) = 0.0_dp
2600 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2601 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2602 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2603 dxpndxk(2, 2) = 0.0_dp
2604 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2605 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2606 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2607 dxpndxk(3, 3) = 0.0_dp
2609 fi = matmul(dsdxpn, dxpndxi)
2610 fj = matmul(dsdxpn, dxpndxj)
2611 fk = matmul(dsdxpn, dxpndxk)
2614 CALL put_derivative(colvar, np + 1, fi)
2615 CALL put_derivative(colvar, np + 2, fj)
2616 CALL put_derivative(colvar, np + 3, fk)
2619 END SUBROUTINE plane_plane_angle_colvar
2629 SUBROUTINE rotation_colvar(colvar, cell, subsys, particles)
2630 TYPE(colvar_type),
POINTER :: colvar
2631 TYPE(cell_type),
POINTER :: cell
2632 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2633 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2634 POINTER :: particles
2637 REAL(
dp) :: a, b, fmod, t0, t1, t2, t3, xdum(3), &
2639 REAL(kind=
dp) :: dp1b1(3), dp1b2(3), dp2b1(3), dp2b2(3), &
2640 ss(3), xp1b1(3), xp1b2(3), xp2b1(3), &
2642 TYPE(particle_list_type),
POINTER :: particles_i
2643 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2645 NULLIFY (particles_i)
2648 IF (
PRESENT(particles))
THEN
2649 my_particles => particles
2651 cpassert(
PRESENT(subsys))
2653 my_particles => particles_i%els
2655 i = colvar%rotation_param%i_at1_bond1
2656 CALL get_coordinates(colvar, i, xp1b1, my_particles)
2657 i = colvar%rotation_param%i_at2_bond1
2658 CALL get_coordinates(colvar, i, xp2b1, my_particles)
2659 i = colvar%rotation_param%i_at1_bond2
2660 CALL get_coordinates(colvar, i, xp1b2, my_particles)
2661 i = colvar%rotation_param%i_at2_bond2
2662 CALL get_coordinates(colvar, i, xp2b2, my_particles)
2664 ss = matmul(cell%h_inv, xp1b1 - xp2b1)
2666 xij = matmul(cell%hmat, ss)
2668 ss = matmul(cell%h_inv, xp1b2 - xp2b2)
2670 xkj = matmul(cell%hmat, ss)
2672 a = sqrt(dot_product(xij, xij))
2673 b = sqrt(dot_product(xkj, xkj))
2675 t1 = 1.0_dp/(a**3.0_dp*b)
2676 t2 = 1.0_dp/(a*b**3.0_dp)
2677 t3 = dot_product(xij, xkj)
2678 colvar%ss = acos(t3*t0)
2679 IF ((abs(colvar%ss) .LT. tolerance_acos) .OR. (abs(colvar%ss -
pi) .LT. tolerance_acos))
THEN
2682 fmod = -1.0_dp/sin(colvar%ss)
2684 dp1b1 = xkj(:)*t0 - xij(:)*t1*t3
2685 dp2b1 = -xkj(:)*t0 + xij(:)*t1*t3
2686 dp1b2 = xij(:)*t0 - xkj(:)*t2*t3
2687 dp2b2 = -xij(:)*t0 + xkj(:)*t2*t3
2690 idum = colvar%rotation_param%i_at1_bond1
2691 CALL put_derivative(colvar, idum, xdum)
2693 idum = colvar%rotation_param%i_at2_bond1
2694 CALL put_derivative(colvar, idum, xdum)
2696 idum = colvar%rotation_param%i_at1_bond2
2697 CALL put_derivative(colvar, idum, xdum)
2699 idum = colvar%rotation_param%i_at2_bond2
2700 CALL put_derivative(colvar, idum, xdum)
2702 END SUBROUTINE rotation_colvar
2713 SUBROUTINE dfunct_colvar(colvar, cell, subsys, particles)
2714 TYPE(colvar_type),
POINTER :: colvar
2715 TYPE(cell_type),
POINTER :: cell
2716 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2717 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2718 POINTER :: particles
2720 INTEGER :: i, j, k, l
2721 REAL(
dp) :: fi(3), fj(3), fk(3), fl(3), r12, r34, &
2722 ss(3), xij(3), xkl(3), xpi(3), xpj(3), &
2724 TYPE(particle_list_type),
POINTER :: particles_i
2725 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2727 NULLIFY (particles_i)
2730 IF (
PRESENT(particles))
THEN
2731 my_particles => particles
2733 cpassert(
PRESENT(subsys))
2735 my_particles => particles_i%els
2737 i = colvar%dfunct_param%i_at_dfunct(1)
2738 j = colvar%dfunct_param%i_at_dfunct(2)
2740 CALL get_coordinates(colvar, i, xpi, my_particles)
2741 CALL get_coordinates(colvar, j, xpj, my_particles)
2742 IF (colvar%dfunct_param%use_pbc)
THEN
2743 ss = matmul(cell%h_inv, xpi - xpj)
2745 xij = matmul(cell%hmat, ss)
2749 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2751 k = colvar%dfunct_param%i_at_dfunct(3)
2752 l = colvar%dfunct_param%i_at_dfunct(4)
2753 CALL get_coordinates(colvar, k, xpk, my_particles)
2754 CALL get_coordinates(colvar, l, xpl, my_particles)
2755 IF (colvar%dfunct_param%use_pbc)
THEN
2756 ss = matmul(cell%h_inv, xpk - xpl)
2758 xkl = matmul(cell%hmat, ss)
2762 r34 = sqrt(xkl(1)**2 + xkl(2)**2 + xkl(3)**2)
2764 colvar%ss = r12 + colvar%dfunct_param%coeff*r34
2767 fk(:) = colvar%dfunct_param%coeff*xkl/r34
2768 fl(:) = -colvar%dfunct_param%coeff*xkl/r34
2769 CALL put_derivative(colvar, 1, fi)
2770 CALL put_derivative(colvar, 2, fj)
2771 CALL put_derivative(colvar, 3, fk)
2772 CALL put_derivative(colvar, 4, fl)
2774 END SUBROUTINE dfunct_colvar
2784 SUBROUTINE angle_colvar(colvar, cell, subsys, particles)
2785 TYPE(colvar_type),
POINTER :: colvar
2786 TYPE(cell_type),
POINTER :: cell
2787 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2788 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2789 POINTER :: particles
2792 REAL(
dp) :: a, b, fi(3), fj(3), fk(3), fmod, ri(3), &
2793 rj(3), rk(3), ss(3), t0, t1, t2, t3, &
2795 TYPE(particle_list_type),
POINTER :: particles_i
2796 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2798 NULLIFY (particles_i)
2801 IF (
PRESENT(particles))
THEN
2802 my_particles => particles
2804 cpassert(
PRESENT(subsys))
2806 my_particles => particles_i%els
2808 i = colvar%angle_param%i_at_angle(1)
2809 j = colvar%angle_param%i_at_angle(2)
2810 k = colvar%angle_param%i_at_angle(3)
2811 CALL get_coordinates(colvar, i, ri, my_particles)
2812 CALL get_coordinates(colvar, j, rj, my_particles)
2813 CALL get_coordinates(colvar, k, rk, my_particles)
2815 ss = matmul(cell%h_inv, ri - rj)
2817 xij = matmul(cell%hmat, ss)
2819 ss = matmul(cell%h_inv, rk - rj)
2821 xkj = matmul(cell%hmat, ss)
2823 a = sqrt(dot_product(xij, xij))
2824 b = sqrt(dot_product(xkj, xkj))
2826 t1 = 1.0_dp/(a**3.0_dp*b)
2827 t2 = 1.0_dp/(a*b**3.0_dp)
2828 t3 = dot_product(xij, xkj)
2829 colvar%ss = acos(t3*t0)
2830 IF ((abs(colvar%ss) .LT. tolerance_acos) .OR. (abs(colvar%ss -
pi) .LT. tolerance_acos))
THEN
2833 fmod = -1.0_dp/sin(colvar%ss)
2835 fi(:) = xkj(:)*t0 - xij(:)*t1*t3
2836 fj(:) = -xkj(:)*t0 + xij(:)*t1*t3 - xij(:)*t0 + xkj(:)*t2*t3
2837 fk(:) = xij(:)*t0 - xkj(:)*t2*t3
2841 CALL put_derivative(colvar, 1, fi)
2842 CALL put_derivative(colvar, 2, fj)
2843 CALL put_derivative(colvar, 3, fk)
2845 END SUBROUTINE angle_colvar
2855 SUBROUTINE dist_colvar(colvar, cell, subsys, particles)
2856 TYPE(colvar_type),
POINTER :: colvar
2857 TYPE(cell_type),
POINTER :: cell
2858 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2859 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2860 POINTER :: particles
2863 REAL(
dp) :: fi(3), fj(3), r12, ss(3), xij(3), &
2865 TYPE(particle_list_type),
POINTER :: particles_i
2866 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2868 NULLIFY (particles_i)
2871 IF (
PRESENT(particles))
THEN
2872 my_particles => particles
2874 cpassert(
PRESENT(subsys))
2876 my_particles => particles_i%els
2878 i = colvar%dist_param%i_at
2879 j = colvar%dist_param%j_at
2880 CALL get_coordinates(colvar, i, xpi, my_particles)
2881 CALL get_coordinates(colvar, j, xpj, my_particles)
2882 ss = matmul(cell%h_inv, xpi - xpj)
2884 xij = matmul(cell%hmat, ss)
2885 SELECT CASE (colvar%dist_param%axis_id)
2904 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2906 IF (colvar%dist_param%sign_d)
THEN
2907 SELECT CASE (colvar%dist_param%axis_id)
2925 CALL put_derivative(colvar, 1, fi)
2926 CALL put_derivative(colvar, 2, fj)
2928 END SUBROUTINE dist_colvar
2939 SUBROUTINE torsion_colvar(colvar, cell, subsys, particles, no_riemann_sheet_op)
2941 TYPE(colvar_type),
POINTER :: colvar
2942 TYPE(cell_type),
POINTER :: cell
2943 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
2944 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
2945 POINTER :: particles
2946 LOGICAL,
INTENT(IN),
OPTIONAL :: no_riemann_sheet_op
2949 LOGICAL :: no_riemann_sheet
2950 REAL(
dp) ::
angle, cosine, dedphi, dedxia, dedxib, dedxic, dedxid, dedxt, dedxu, dedyia, &
2951 dedyib, dedyic, dedyid, dedyt, dedyu, dedzia, dedzib, dedzic, dedzid, dedzt, dedzu, dt, &
2952 e, ftmp(3), o0, rcb, rt2, rtmp(3), rtru, ru2, sine, ss(3), xba, xca, xcb, xdb, xdc, xt, &
2953 xtu, xu, yba, yca, ycb, ydb, ydc, yt, ytu, yu, zba, zca, zcb, zdb, zdc, zt, ztu, zu
2954 REAL(
dp),
DIMENSION(3, 4) :: rr
2955 TYPE(particle_list_type),
POINTER :: particles_i
2956 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
2958 NULLIFY (particles_i)
2960 IF (
PRESENT(particles))
THEN
2961 my_particles => particles
2963 cpassert(
PRESENT(subsys))
2965 my_particles => particles_i%els
2967 no_riemann_sheet = .false.
2968 IF (
PRESENT(no_riemann_sheet_op)) no_riemann_sheet = no_riemann_sheet_op
2970 i = colvar%torsion_param%i_at_tors(ii)
2971 CALL get_coordinates(colvar, i, rtmp, my_particles)
2972 rr(:, ii) = rtmp(1:3)
2974 o0 = colvar%torsion_param%o0
2976 ss = matmul(cell%h_inv, rr(:, 2) - rr(:, 1))
2978 ss = matmul(cell%hmat, ss)
2983 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 2))
2985 ss = matmul(cell%hmat, ss)
2990 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 3))
2992 ss = matmul(cell%hmat, ss)
2997 xt = yba*zcb - ycb*zba
2998 yt = zba*xcb - zcb*xba
2999 zt = xba*ycb - xcb*yba
3000 xu = ycb*zdc - ydc*zcb
3001 yu = zcb*xdc - zdc*xcb
3002 zu = xcb*ydc - xdc*ycb
3006 rt2 = xt*xt + yt*yt + zt*zt
3007 ru2 = xu*xu + yu*yu + zu*zu
3008 rtru = sqrt(rt2*ru2)
3009 IF (rtru .NE. 0.0_dp)
THEN
3010 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb)
3011 cosine = (xt*xu + yt*yu + zt*zu)/rtru
3012 sine = (xcb*xtu + ycb*ytu + zcb*ztu)/(rcb*rtru)
3013 cosine = min(1.0_dp, max(-1.0_dp, cosine))
3014 angle = acos(cosine)
3018 dt = mod(2.0e4_dp*
pi + dt - o0, 2.0_dp*
pi)
3019 IF (dt .GT.
pi) dt = dt - 2.0_dp*
pi
3021 colvar%torsion_param%o0 = dt
3031 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 1))
3033 ss = matmul(cell%hmat, ss)
3038 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 2))
3040 ss = matmul(cell%hmat, ss)
3045 dedxt = dedphi*(yt*zcb - ycb*zt)/(rt2*rcb)
3046 dedyt = dedphi*(zt*xcb - zcb*xt)/(rt2*rcb)
3047 dedzt = dedphi*(xt*ycb - xcb*yt)/(rt2*rcb)
3048 dedxu = -dedphi*(yu*zcb - ycb*zu)/(ru2*rcb)
3049 dedyu = -dedphi*(zu*xcb - zcb*xu)/(ru2*rcb)
3050 dedzu = -dedphi*(xu*ycb - xcb*yu)/(ru2*rcb)
3054 dedxia = zcb*dedyt - ycb*dedzt
3055 dedyia = xcb*dedzt - zcb*dedxt
3056 dedzia = ycb*dedxt - xcb*dedyt
3057 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu
3058 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu
3059 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu
3060 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu
3061 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu
3062 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu
3063 dedxid = zcb*dedyu - ycb*dedzu
3064 dedyid = xcb*dedzu - zcb*dedxu
3065 dedzid = ycb*dedxu - xcb*dedyu
3082 IF (no_riemann_sheet) colvar%ss = atan2(sin(e), cos(e))
3086 CALL put_derivative(colvar, 1, ftmp)
3090 CALL put_derivative(colvar, 2, ftmp)
3094 CALL put_derivative(colvar, 3, ftmp)
3098 CALL put_derivative(colvar, 4, ftmp)
3099 END SUBROUTINE torsion_colvar
3108 SUBROUTINE qparm_colvar(colvar, cell, subsys, particles)
3109 TYPE(colvar_type),
POINTER :: colvar
3110 TYPE(cell_type),
POINTER :: cell
3111 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
3112 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
3113 POINTER :: particles
3115 INTEGER :: aa, bb, cc, i, idim, ii, j, jj, l, mm, &
3116 n_atoms_from, n_atoms_to, ncells(3)
3117 LOGICAL :: include_images
3118 REAL(kind=
dp) :: denominator_tolerance, fact, ftmp(3), im_qlm, inv_n_atoms_from, nbond, &
3119 pre_fac, ql, qparm, r1cut, rcut, re_qlm, rij, rij_shift, shift(3), ss(3), ss0(3), xij(3), &
3121 REAL(kind=
dp),
DIMENSION(3) :: d_im_qlm_dxi, d_nbond_dxi, d_ql_dxi, &
3122 d_re_qlm_dxi, xpi, xpj
3123 TYPE(particle_list_type),
POINTER :: particles_i
3124 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
3130 n_atoms_to = colvar%qparm_param%n_atoms_to
3131 n_atoms_from = colvar%qparm_param%n_atoms_from
3132 rcut = colvar%qparm_param%rcut
3133 l = colvar%qparm_param%l
3134 r1cut = colvar%qparm_param%rstart
3135 include_images = colvar%qparm_param%include_images
3136 NULLIFY (particles_i)
3138 IF (
PRESENT(particles))
THEN
3139 my_particles => particles
3141 cpassert(
PRESENT(subsys))
3143 my_particles => particles_i%els
3145 cpassert(r1cut .LT. rcut)
3146 denominator_tolerance = 1.0e-8_dp
3153 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
3154 DO ii = 1, n_atoms_from
3155 i = colvar%qparm_param%i_at_from(ii)
3156 CALL get_coordinates(colvar, i, xpi, my_particles)
3159 d_ql_dxi(:) = 0.0_dp
3165 d_re_qlm_dxi(:) = 0.0_dp
3166 d_im_qlm_dxi(:) = 0.0_dp
3167 d_nbond_dxi(:) = 0.0_dp
3169 jloop:
DO jj = 1, n_atoms_to
3171 j = colvar%qparm_param%i_at_to(jj)
3172 CALL get_coordinates(colvar, j, xpj, my_particles)
3174 IF (include_images)
THEN
3176 cpassert(cell%orthorhombic)
3180 xij(:) = xpj(:) - xpi(:)
3181 ss = matmul(cell%h_inv, xij)
3187 shift(idim) = 1.0_dp
3188 xij_shift = matmul(cell%hmat, shift)
3189 rij_shift = sqrt(dot_product(xij_shift, xij_shift))
3190 ncells(idim) = floor(rcut/rij_shift - 0.5)
3195 DO aa = -ncells(1), ncells(1)
3196 DO bb = -ncells(2), ncells(2)
3197 DO cc = -ncells(3), ncells(3)
3199 IF (i == j .AND. aa .EQ. 0 .AND. bb .EQ. 0 .AND. cc .EQ. 0) cycle
3200 shift(1) = real(aa, kind=
dp)
3201 shift(2) = real(bb, kind=
dp)
3202 shift(3) = real(cc, kind=
dp)
3203 xij = matmul(cell%hmat, ss0(:) + shift(:))
3204 rij = sqrt(dot_product(xij, xij))
3210 IF (rij > rcut) cycle
3213 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3214 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3215 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3223 IF (i == j) cycle jloop
3224 xij(:) = xpj(:) - xpi(:)
3225 rij = sqrt(dot_product(xij, xij))
3226 IF (rij > rcut) cycle
3229 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3230 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3231 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3246 IF (nbond .LT. denominator_tolerance)
THEN
3247 cpwarn(
"QPARM: number of neighbors is very close to zero!")
3250 d_nbond_dxi(:) = d_nbond_dxi(:)/nbond
3251 re_qlm = re_qlm/nbond
3252 d_re_qlm_dxi(:) = d_re_qlm_dxi(:)/nbond - d_nbond_dxi(:)*re_qlm
3253 im_qlm = im_qlm/nbond
3254 d_im_qlm_dxi(:) = d_im_qlm_dxi(:)/nbond - d_nbond_dxi(:)*im_qlm
3256 ql = ql + fact*(re_qlm*re_qlm + im_qlm*im_qlm)
3257 d_ql_dxi(:) = d_ql_dxi(:) &
3258 + fact*2.0_dp*(re_qlm*d_re_qlm_dxi(:) + im_qlm*d_im_qlm_dxi(:))
3262 pre_fac = (4.0_dp*
pi)/(2.0_dp*l + 1)
3264 qparm = qparm + sqrt(pre_fac*ql)
3265 ftmp(:) = 0.5_dp*sqrt(pre_fac/ql)*d_ql_dxi(:)
3267 ftmp(:) = -1.0_dp*ftmp(:)
3269 CALL put_derivative(colvar, ii, ftmp)
3273 colvar%ss = qparm*inv_n_atoms_from
3274 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
3280 END SUBROUTINE qparm_colvar
3298 SUBROUTINE accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3299 denominator_tolerance, ll, mm, nbond, re_qlm, im_qlm, &
3300 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3302 REAL(kind=
dp),
INTENT(IN) :: xij(3), rij, rcut, r1cut, &
3303 denominator_tolerance
3304 INTEGER,
INTENT(IN) :: ll, mm
3305 REAL(kind=
dp),
INTENT(INOUT) :: nbond, re_qlm, im_qlm, d_re_qlm_dxi(3), &
3306 d_im_qlm_dxi(3), d_nbond_dxi(3)
3308 REAL(kind=
dp) :: bond, costheta, dplm, dylm, exp0, &
3309 exp_fac, fi, plm, pre_fac, sqrt_c1
3310 REAL(kind=
dp),
DIMENSION(3) :: dcostheta, dfi
3315 IF (rij .GT. rcut)
THEN
3320 IF (rij .LT. r1cut)
THEN
3324 exp0 = exp((r1cut - rcut)/(rij - rcut) - (r1cut - rcut)/(r1cut - rij))
3325 bond = 1.0_dp/(1.0_dp + exp0)
3326 exp_fac = ((rcut - r1cut)/(rij - rcut)**2 + (rcut - r1cut)/(r1cut - rij)**2)*exp0/(1.0_dp + exp0)**2
3329 IF (bond > 1.0_dp)
THEN
3330 cpabort(
"bond > 1.0_dp")
3333 nbond = nbond + bond
3334 IF (abs(xij(1)) .LT. denominator_tolerance &
3335 .AND. abs(xij(2)) .LT. denominator_tolerance)
THEN
3338 fi = atan2(xij(2), xij(1))
3341 costheta = xij(3)/rij
3342 IF (costheta > 1.0_dp) costheta = 1.0_dp
3343 IF (costheta < -1.0_dp) costheta = -1.0_dp
3348 IF ((ll + abs(mm)) >
maxfac)
THEN
3349 cpabort(
"(l+m) > maxfac")
3352 sqrt_c1 = sqrt(((2*ll + 1)*
fac(ll - abs(mm)))/(4*
pi*
fac(ll + abs(mm))))
3353 pre_fac = bond*sqrt_c1
3361 re_qlm = re_qlm + pre_fac*plm*cos(mm*fi)
3362 im_qlm = im_qlm + pre_fac*plm*sin(mm*fi)
3367 dcostheta(:) = xij(:)*xij(3)/(rij**3)
3368 dcostheta(3) = dcostheta(3) - 1.0_dp/rij
3372 dfi(1) = xij(2)/(xij(1)**2 + xij(2)**2)
3373 dfi(2) = -xij(1)/(xij(1)**2 + xij(2)**2)
3375 d_re_qlm_dxi(:) = d_re_qlm_dxi(:) &
3376 + exp_fac*sqrt_c1*plm*cos(mm*fi)*xij(:)/rij &
3377 + dylm*dcostheta(:)*cos(mm*fi) &
3378 + pre_fac*plm*mm*(-1.0_dp)*sin(mm*fi)*dfi(:)
3379 d_im_qlm_dxi(:) = d_im_qlm_dxi(:) &
3380 + exp_fac*sqrt_c1*plm*sin(mm*fi)*xij(:)/rij &
3381 + dylm*dcostheta(:)*sin(mm*fi) &
3382 + pre_fac*plm*mm*(+1.0_dp)*cos(mm*fi)*dfi(:)
3383 d_nbond_dxi(:) = d_nbond_dxi(:) + exp_fac*xij(:)/rij
3385 END SUBROUTINE accumulate_qlm_over_neigbors
3397 SUBROUTINE hydronium_shell_colvar(colvar, cell, subsys, particles)
3398 TYPE(colvar_type),
POINTER :: colvar
3399 TYPE(cell_type),
POINTER :: cell
3400 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
3401 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
3402 POINTER :: particles
3404 INTEGER :: i, ii, j, jj, n_hydrogens, n_oxygens, &
3405 pm, poh, poo, qm, qoh, qoo
3406 REAL(
dp) :: drji, fscalar, invden, lambda, nh, num, &
3407 qtot, rji(3), roh, roo, rrel
3408 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: m, noh, noo, qloc
3409 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dm, dnoh, dnoo
3410 REAL(
dp),
DIMENSION(3) :: rpi, rpj
3411 TYPE(particle_list_type),
POINTER :: particles_i
3412 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
3414 n_oxygens = colvar%hydronium_shell_param%n_oxygens
3415 n_hydrogens = colvar%hydronium_shell_param%n_hydrogens
3416 nh = colvar%hydronium_shell_param%nh
3417 poh = colvar%hydronium_shell_param%poh
3418 qoh = colvar%hydronium_shell_param%qoh
3419 poo = colvar%hydronium_shell_param%poo
3420 qoo = colvar%hydronium_shell_param%qoo
3421 roo = colvar%hydronium_shell_param%roo
3422 roh = colvar%hydronium_shell_param%roh
3423 lambda = colvar%hydronium_shell_param%lambda
3424 pm = colvar%hydronium_shell_param%pm
3425 qm = colvar%hydronium_shell_param%qm
3427 NULLIFY (particles_i)
3429 IF (
PRESENT(particles))
THEN
3430 my_particles => particles
3432 cpassert(
PRESENT(subsys))
3434 my_particles => particles_i%els
3437 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3438 ALLOCATE (noh(n_oxygens))
3439 ALLOCATE (m(n_oxygens))
3440 ALLOCATE (dm(3, n_hydrogens, n_oxygens))
3442 ALLOCATE (dnoo(3, n_oxygens, n_oxygens))
3443 ALLOCATE (noo(n_oxygens))
3445 ALLOCATE (qloc(n_oxygens))
3455 DO ii = 1, n_oxygens
3456 i = colvar%hydronium_shell_param%i_oxygens(ii)
3457 rpi(:) = my_particles(i)%r(1:3)
3459 DO jj = 1, n_hydrogens
3460 j = colvar%hydronium_shell_param%i_hydrogens(jj)
3461 rpj(:) = my_particles(j)%r(1:3)
3462 rji =
pbc(rpj, rpi, cell)
3463 drji = sqrt(sum(rji**2))
3465 num = (1.0_dp - rrel**poh)
3466 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3467 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3468 noh(ii) = noh(ii) + num*invden
3469 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3470 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3471 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3474 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3475 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3476 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3479 m(ii) = 1.0_dp - (1.0_dp - (noh(ii)/nh)**pm)/ &
3480 (1.0_dp - (noh(ii)/nh)**qm)
3483 DO jj = 1, n_oxygens
3485 j = colvar%hydronium_shell_param%i_oxygens(jj)
3486 rpj(:) = my_particles(j)%r(1:3)
3487 rji =
pbc(rpj, rpi, cell)
3488 drji = sqrt(sum(rji**2))
3490 num = (1.0_dp - rrel**poo)
3491 invden = 1.0_dp/(1.0_dp - rrel**qoo)
3492 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3493 noo(ii) = noo(ii) + num*invden
3494 fscalar = ((-poo*(rrel**(poo - 1))*invden) &
3495 + num*(invden)**2*qoo*(rrel**(qoo - 1)))/(drji*roo)
3496 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3499 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
3500 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drji)
3501 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3508 DO ii = 1, n_oxygens
3509 qloc(ii) = exp(lambda*m(ii)*noo(ii))
3510 qtot = qtot + qloc(ii)
3513 DO ii = 1, n_oxygens
3515 DO jj = 1, n_hydrogens
3516 dm(1:3, jj, ii) = (pm*((noh(ii)/nh)**(pm - 1))*dnoh(1:3, jj, ii))/nh/ &
3517 (1.0_dp - (noh(ii)/nh)**qm) - &
3518 (1.0_dp - (noh(ii)/nh)**pm)/ &
3519 ((1.0_dp - (noh(ii)/nh)**qm)**2)* &
3520 qm*dnoh(1:3, jj, ii)*(noh(ii)/nh)**(qm - 1)/nh
3522 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3523 colvar%dsdr(1:3, n_oxygens + jj) = colvar%dsdr(1:3, n_oxygens + jj) &
3524 - qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3527 DO jj = 1, n_oxygens
3528 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3529 colvar%dsdr(1:3, jj) = colvar%dsdr(1:3, jj) &
3530 - qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3534 colvar%ss = log(qtot)/lambda
3543 END SUBROUTINE hydronium_shell_colvar
3554 SUBROUTINE hydronium_dist_colvar(colvar, cell, subsys, particles)
3555 TYPE(colvar_type),
POINTER :: colvar
3556 TYPE(cell_type),
POINTER :: cell
3557 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
3558 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
3559 POINTER :: particles
3561 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3562 n_oxygens, offseth, pf, pm, poh, qf, &
3564 REAL(
dp) :: drji, drki, fscalar, invden, lambda, nh, nn, num, rion, rion_den, rion_num, &
3565 rji(3), rki(3), roh, rrel, sum_expfac_f, sum_expfac_noh
3566 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac_f, dexpfac_noh, df, dm, &
3567 expfac_f, expfac_f_rki, expfac_noh, f, &
3569 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_f_rki
3570 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rki, dnoh
3571 REAL(
dp),
DIMENSION(3) :: rpi, rpj, rpk
3572 TYPE(particle_list_type),
POINTER :: particles_i
3573 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
3575 n_oxygens = colvar%hydronium_dist_param%n_oxygens
3576 n_hydrogens = colvar%hydronium_dist_param%n_hydrogens
3577 poh = colvar%hydronium_dist_param%poh
3578 qoh = colvar%hydronium_dist_param%qoh
3579 roh = colvar%hydronium_dist_param%roh
3580 pm = colvar%hydronium_dist_param%pm
3581 qm = colvar%hydronium_dist_param%qm
3582 nh = colvar%hydronium_dist_param%nh
3583 pf = colvar%hydronium_dist_param%pf
3584 qf = colvar%hydronium_dist_param%qf
3585 nn = colvar%hydronium_dist_param%nn
3586 lambda = colvar%hydronium_dist_param%lambda
3588 NULLIFY (particles_i)
3590 IF (
PRESENT(particles))
THEN
3591 my_particles => particles
3593 cpassert(
PRESENT(subsys))
3595 my_particles => particles_i%els
3598 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3599 ALLOCATE (noh(n_oxygens))
3600 ALLOCATE (m(n_oxygens), dm(n_oxygens))
3601 ALLOCATE (f(n_oxygens), df(n_oxygens))
3602 ALLOCATE (expfac_noh(n_oxygens), dexpfac_noh(n_oxygens))
3603 ALLOCATE (expfac_f(n_oxygens), dexpfac_f(n_oxygens))
3604 ALLOCATE (ddist_rki(3, n_oxygens, n_oxygens))
3605 ALLOCATE (expfac_f_rki(n_oxygens))
3606 ALLOCATE (dexpfac_f_rki(n_oxygens, n_oxygens))
3618 sum_expfac_noh = 0._dp
3619 sum_expfac_f = 0._dp
3621 expfac_f_rki = 0._dp
3622 dexpfac_f_rki = 0._dp
3625 DO ii = 1, n_oxygens
3626 i = colvar%hydronium_dist_param%i_oxygens(ii)
3627 rpi(:) = my_particles(i)%r(1:3)
3628 DO jj = 1, n_hydrogens
3629 j = colvar%hydronium_dist_param%i_hydrogens(jj)
3630 rpj(:) = my_particles(j)%r(1:3)
3631 rji =
pbc(rpj, rpi, cell)
3632 drji = sqrt(sum(rji**2))
3634 num = (1.0_dp - rrel**poh)
3635 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3636 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3637 noh(ii) = noh(ii) + num*invden
3638 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3639 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3640 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3643 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3644 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3645 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3651 DO ii = 1, n_oxygens
3652 num = 1.0_dp - (noh(ii)/nh)**pm
3653 invden = 1.0_dp/(1.0_dp - (noh(ii)/nh)**qm)
3654 m(ii) = 1.0_dp - num*invden
3655 dm(ii) = (pm*(noh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
3656 (noh(ii)/nh)**(qm - 1))/nh
3657 expfac_noh(ii) = exp(lambda*noh(ii))
3658 dexpfac_noh(ii) = lambda*expfac_noh(ii)
3659 sum_expfac_noh = sum_expfac_noh + expfac_noh(ii)
3663 DO ii = 1, n_oxygens
3664 i = colvar%hydronium_dist_param%i_oxygens(ii)
3665 num = 1.0_dp - (noh(ii)/nn)**pf
3666 invden = 1.0_dp/(1.0_dp - (noh(ii)/nn)**qf)
3668 df(ii) = (-pf*(noh(ii)/nn)**(pf - 1)*invden + qf*num*(invden**2)* &
3669 (noh(ii)/nn)**(qf - 1))/nn
3670 expfac_f(ii) = exp(lambda*f(ii))
3671 dexpfac_f(ii) = lambda*expfac_f(ii)
3672 sum_expfac_f = sum_expfac_f + expfac_f(ii)
3676 DO ii = 1, n_oxygens
3677 i = colvar%hydronium_dist_param%i_oxygens(ii)
3678 rpi(:) = my_particles(i)%r(1:3)
3679 DO kk = 1, n_oxygens
3681 k = colvar%hydronium_dist_param%i_oxygens(kk)
3682 rpk(:) = my_particles(k)%r(1:3)
3683 rki =
pbc(rpk, rpi, cell)
3684 drki = sqrt(sum(rki**2))
3685 expfac_f_rki(ii) = expfac_f_rki(ii) + drki*expfac_f(kk)
3686 ddist_rki(1:3, kk, ii) = rki(1:3)/drki
3687 dexpfac_f_rki(kk, ii) = drki*dexpfac_f(kk)
3689 rion_num = rion_num + m(ii)*expfac_noh(ii)*expfac_f_rki(ii)
3693 rion_den = sum_expfac_noh*sum_expfac_f
3694 rion = rion_num/rion_den
3699 DO ii = 1, n_oxygens
3700 DO jj = 1, n_hydrogens
3701 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3702 + dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3703 *expfac_f_rki(ii)/rion_den
3704 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3705 - dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3706 *expfac_f_rki(ii)/rion_den
3707 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3708 + m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3709 *expfac_f_rki(ii)/rion_den
3710 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3711 - m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3712 *expfac_f_rki(ii)/rion_den
3714 DO kk = 1, n_oxygens
3716 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3717 - m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3718 *expfac_f(kk)/rion_den
3719 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3720 + m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3721 *expfac_f(kk)/rion_den
3722 DO jj = 1, n_hydrogens
3723 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3724 + m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3725 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3726 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3727 - m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3728 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3733 DO ii = 1, n_oxygens
3734 DO jj = 1, n_hydrogens
3735 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3736 - rion_num*sum_expfac_f*dexpfac_noh(ii) &
3737 *dnoh(1:3, jj, ii)/(rion_den**2)
3738 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3739 + rion_num*sum_expfac_f*dexpfac_noh(ii) &
3740 *dnoh(1:3, jj, ii)/(rion_den**2)
3741 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3742 - rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3743 *dnoh(1:3, jj, ii)/(rion_den**2)
3744 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3745 + rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3746 *dnoh(1:3, jj, ii)/(rion_den**2)
3750 DEALLOCATE (noh, m, f, expfac_noh, expfac_f)
3751 DEALLOCATE (dnoh, dm, df, dexpfac_noh, dexpfac_f)
3752 DEALLOCATE (ddist_rki, expfac_f_rki, dexpfac_f_rki)
3754 END SUBROUTINE hydronium_dist_colvar
3767 SUBROUTINE acid_hyd_dist_colvar(colvar, cell, subsys, particles)
3768 TYPE(colvar_type),
POINTER :: colvar
3769 TYPE(cell_type),
POINTER :: cell
3770 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
3771 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
3772 POINTER :: particles
3774 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3775 n_oxygens_acid, n_oxygens_water, &
3776 offseth, offseto, paoh, pcut, pwoh, &
3778 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac, expfac, nwoh
3779 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_rik
3780 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rik, dnaoh, dnwoh
3781 REAL(kind=
dp) :: dfcut, drik, drji, drjk, fbrace, fcut, fscalar, invden, invden_cut, lambda, &
3782 naoh, nc, num, num_cut, raoh, rik(3), rion, rion_den, rion_num, rji(3), rjk(3), rpi(3), &
3783 rpj(3), rpk(3), rrel, rwoh
3784 TYPE(particle_list_type),
POINTER :: particles_i
3785 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
3787 NULLIFY (my_particles, particles_i)
3789 n_oxygens_water = colvar%acid_hyd_dist_param%n_oxygens_water
3790 n_oxygens_acid = colvar%acid_hyd_dist_param%n_oxygens_acid
3791 n_hydrogens = colvar%acid_hyd_dist_param%n_hydrogens
3792 pwoh = colvar%acid_hyd_dist_param%pwoh
3793 qwoh = colvar%acid_hyd_dist_param%qwoh
3794 paoh = colvar%acid_hyd_dist_param%paoh
3795 qaoh = colvar%acid_hyd_dist_param%qaoh
3796 pcut = colvar%acid_hyd_dist_param%pcut
3797 qcut = colvar%acid_hyd_dist_param%qcut
3798 rwoh = colvar%acid_hyd_dist_param%rwoh
3799 raoh = colvar%acid_hyd_dist_param%raoh
3800 nc = colvar%acid_hyd_dist_param%nc
3801 lambda = colvar%acid_hyd_dist_param%lambda
3802 ALLOCATE (expfac(n_oxygens_water))
3803 ALLOCATE (nwoh(n_oxygens_water))
3804 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
3805 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
3806 ALLOCATE (dexpfac(n_oxygens_water))
3807 ALLOCATE (ddist_rik(3, n_oxygens_water, n_oxygens_acid))
3808 ALLOCATE (dexpfac_rik(n_oxygens_water, n_oxygens_acid))
3813 dnaoh(:, :, :) = 0._dp
3814 dnwoh(:, :, :) = 0._dp
3815 ddist_rik(:, :, :) = 0._dp
3817 dexpfac_rik(:, :) = 0._dp
3820 IF (
PRESENT(particles))
THEN
3821 my_particles => particles
3823 cpassert(
PRESENT(subsys))
3825 my_particles => particles_i%els
3829 DO ii = 1, n_oxygens_water
3830 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3831 rpi(:) = my_particles(i)%r(1:3)
3832 DO jj = 1, n_hydrogens
3833 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3834 rpj(:) = my_particles(j)%r(1:3)
3835 rji =
pbc(rpj, rpi, cell)
3836 drji = sqrt(sum(rji**2))
3838 num = 1.0_dp - rrel**pwoh
3839 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
3840 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3841 nwoh(ii) = nwoh(ii) + num*invden
3842 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
3843 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
3844 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3847 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
3848 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
3849 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3852 expfac(ii) = exp(lambda*nwoh(ii))
3853 dexpfac(ii) = lambda*expfac(ii)
3854 rion_den = rion_den + expfac(ii)
3858 DO kk = 1, n_oxygens_acid
3859 k = colvar%acid_hyd_dist_param%i_oxygens_acid(kk)
3860 rpk(:) = my_particles(k)%r(1:3)
3861 DO ii = 1, n_oxygens_water
3862 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3863 rpi(:) = my_particles(i)%r(1:3)
3864 rik =
pbc(rpi, rpk, cell)
3865 drik = sqrt(sum(rik**2))
3866 rion_num = rion_num + drik*expfac(ii)
3867 ddist_rik(1:3, ii, kk) = rik(1:3)/drik
3868 dexpfac_rik(ii, kk) = drik*dexpfac(ii)
3873 DO kk = 1, n_oxygens_acid
3874 k = colvar%acid_hyd_dist_param%i_oxygens_acid(kk)
3875 rpk(:) = my_particles(k)%r(1:3)
3876 DO jj = 1, n_hydrogens
3877 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3878 rpj(:) = my_particles(j)%r(1:3)
3879 rjk =
pbc(rpj, rpk, cell)
3880 drjk = sqrt(sum(rjk**2))
3882 num = 1.0_dp - rrel**paoh
3883 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
3884 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3885 naoh = naoh + num*invden
3886 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
3887 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
3888 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3891 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
3892 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
3893 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3897 num_cut = 1.0_dp - (naoh/nc)**pcut
3898 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
3899 fcut = num_cut*invden_cut
3903 fbrace = rion_num/rion_den/n_oxygens_acid
3908 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
3909 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
3910 offseto = n_oxygens_water
3911 offseth = n_oxygens_water + n_oxygens_acid
3912 DO kk = 1, n_oxygens_acid
3913 DO jj = 1, n_hydrogens
3914 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3915 + dfcut*dnaoh(1:3, jj, kk)*fbrace
3916 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3917 - dfcut*dnaoh(1:3, jj, kk)*fbrace
3923 DO kk = 1, n_oxygens_acid
3924 DO ii = 1, n_oxygens_water
3925 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3926 + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3928 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3929 - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3931 DO jj = 1, n_hydrogens
3932 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3933 + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3935 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3936 - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3942 DO ii = 1, n_oxygens_water
3943 DO jj = 1, n_hydrogens
3944 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3945 - fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3946 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3947 + fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3951 END SUBROUTINE acid_hyd_dist_colvar
3964 SUBROUTINE acid_hyd_shell_colvar(colvar, cell, subsys, particles)
3965 TYPE(colvar_type),
POINTER :: colvar
3966 TYPE(cell_type),
POINTER :: cell
3967 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
3968 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
3969 POINTER :: particles
3971 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, n_oxygens_acid, n_oxygens_water, offseth, &
3972 offseto, paoh, pcut, pm, poo, pwoh, qaoh, qcut, qm, qoo, qwoh, tt
3973 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dm, m, noo, nwoh, qloc
3974 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dnaoh, dnoo, dnwoh
3975 REAL(kind=
dp) :: dfcut, drji, drjk, drki, fcut, fscalar, invden, invden_cut, lambda, naoh, &
3976 nc, nh, num, num_cut, qsol, qtot, raoh, rji(3), rjk(3), rki(3), roo, rpi(3), rpj(3), &
3978 TYPE(particle_list_type),
POINTER :: particles_i
3979 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
3981 NULLIFY (my_particles, particles_i)
3983 n_oxygens_water = colvar%acid_hyd_shell_param%n_oxygens_water
3984 n_oxygens_acid = colvar%acid_hyd_shell_param%n_oxygens_acid
3985 n_hydrogens = colvar%acid_hyd_shell_param%n_hydrogens
3986 pwoh = colvar%acid_hyd_shell_param%pwoh
3987 qwoh = colvar%acid_hyd_shell_param%qwoh
3988 paoh = colvar%acid_hyd_shell_param%paoh
3989 qaoh = colvar%acid_hyd_shell_param%qaoh
3990 poo = colvar%acid_hyd_shell_param%poo
3991 qoo = colvar%acid_hyd_shell_param%qoo
3992 pm = colvar%acid_hyd_shell_param%pm
3993 qm = colvar%acid_hyd_shell_param%qm
3994 pcut = colvar%acid_hyd_shell_param%pcut
3995 qcut = colvar%acid_hyd_shell_param%qcut
3996 rwoh = colvar%acid_hyd_shell_param%rwoh
3997 raoh = colvar%acid_hyd_shell_param%raoh
3998 roo = colvar%acid_hyd_shell_param%roo
3999 nc = colvar%acid_hyd_shell_param%nc
4000 nh = colvar%acid_hyd_shell_param%nh
4001 lambda = colvar%acid_hyd_shell_param%lambda
4002 ALLOCATE (nwoh(n_oxygens_water))
4003 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
4004 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
4005 ALLOCATE (m(n_oxygens_water))
4006 ALLOCATE (dm(n_oxygens_water))
4007 ALLOCATE (noo(n_oxygens_water))
4008 ALLOCATE (dnoo(3, n_oxygens_water + n_oxygens_acid, n_oxygens_water))
4009 ALLOCATE (qloc(n_oxygens_water))
4013 dnaoh(:, :, :) = 0._dp
4014 dnwoh(:, :, :) = 0._dp
4015 dnoo(:, :, :) = 0._dp
4021 IF (
PRESENT(particles))
THEN
4022 my_particles => particles
4024 cpassert(
PRESENT(subsys))
4026 my_particles => particles_i%els
4030 DO ii = 1, n_oxygens_water
4031 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4032 rpi(:) = my_particles(i)%r(1:3)
4033 DO jj = 1, n_hydrogens
4034 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4035 rpj(:) = my_particles(j)%r(1:3)
4036 rji =
pbc(rpj, rpi, cell)
4037 drji = sqrt(sum(rji**2))
4039 num = 1.0_dp - rrel**pwoh
4040 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
4041 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4042 nwoh(ii) = nwoh(ii) + num*invden
4043 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
4044 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
4045 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4048 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
4049 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
4050 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4056 DO ii = 1, n_oxygens_water
4057 num = 1.0_dp - (nwoh(ii)/nh)**pm
4058 invden = 1.0_dp/(1.0_dp - (nwoh(ii)/nh)**qm)
4059 m(ii) = 1.0_dp - num*invden
4060 dm(ii) = (pm*(nwoh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
4061 (nwoh(ii)/nh)**(qm - 1))/nh
4065 DO ii = 1, n_oxygens_water
4066 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4067 rpi(:) = my_particles(i)%r(1:3)
4068 DO kk = 1, n_oxygens_water + n_oxygens_acid
4070 IF (kk <= n_oxygens_water)
THEN
4071 k = colvar%acid_hyd_shell_param%i_oxygens_water(kk)
4072 rpk(:) = my_particles(k)%r(1:3)
4074 tt = kk - n_oxygens_water
4075 k = colvar%acid_hyd_shell_param%i_oxygens_acid(tt)
4076 rpk(:) = my_particles(k)%r(1:3)
4078 rki =
pbc(rpk, rpi, cell)
4079 drki = sqrt(sum(rki**2))
4081 num = 1.0_dp - rrel**poo
4082 invden = 1.0_dp/(1.0_dp - rrel**qoo)
4083 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4084 noo(ii) = noo(ii) + num*invden
4085 fscalar = (-poo*(rrel**(poo - 1))*invden &
4086 + num*(invden**2)*qoo*(rrel**(qoo - 1)))/(drki*roo)
4087 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4090 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
4091 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drki)
4092 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4098 DO kk = 1, n_oxygens_acid
4099 k = colvar%acid_hyd_shell_param%i_oxygens_acid(kk)
4100 rpk(:) = my_particles(k)%r(1:3)
4101 DO jj = 1, n_hydrogens
4102 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4103 rpj(:) = my_particles(j)%r(1:3)
4104 rjk =
pbc(rpj, rpk, cell)
4105 drjk = sqrt(sum(rjk**2))
4107 num = 1.0_dp - rrel**paoh
4108 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
4109 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4110 naoh = naoh + num*invden
4111 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
4112 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
4113 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4116 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
4117 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
4118 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4122 num_cut = 1.0_dp - (naoh/nc)**pcut
4123 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
4124 fcut = num_cut*invden_cut
4127 DO ii = 1, n_oxygens_water
4128 qloc(ii) = exp(lambda*m(ii)*noo(ii))
4129 qtot = qtot + qloc(ii)
4131 qsol = log(qtot)/lambda
4132 colvar%ss = fcut*qsol
4135 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
4136 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
4137 offseto = n_oxygens_water
4138 offseth = n_oxygens_water + n_oxygens_acid
4139 DO kk = 1, n_oxygens_acid
4140 DO jj = 1, n_hydrogens
4141 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
4142 + dfcut*dnaoh(1:3, jj, kk)*qsol
4143 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4144 - dfcut*dnaoh(1:3, jj, kk)*qsol
4150 DO ii = 1, n_oxygens_water
4151 fscalar = fcut*qloc(ii)*dm(ii)*noo(ii)/qtot
4152 DO jj = 1, n_hydrogens
4153 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
4154 + fscalar*dnwoh(1:3, jj, ii)
4155 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4156 - fscalar*dnwoh(1:3, jj, ii)
4160 DO ii = 1, n_oxygens_water
4161 fscalar = fcut*qloc(ii)*m(ii)/qtot
4162 DO kk = 1, n_oxygens_water + n_oxygens_acid
4164 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + fscalar*dnoo(1:3, kk, ii)
4165 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) - fscalar*dnoo(1:3, kk, ii)
4169 END SUBROUTINE acid_hyd_shell_colvar
4181 SUBROUTINE coord_colvar(colvar, cell, subsys, particles)
4182 TYPE(colvar_type),
POINTER :: colvar
4183 TYPE(cell_type),
POINTER :: cell
4184 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
4185 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
4186 POINTER :: particles
4188 INTEGER :: i, ii, j, jj, k, kk, n_atoms_from, &
4189 n_atoms_to_a, n_atoms_to_b, p_a, p_b, &
4191 REAL(
dp) :: dfunc_ij, dfunc_jk, func_ij, func_jk, func_k, inv_n_atoms_from, invden_ij, &
4192 invden_jk, ncoord, num_ij, num_jk, r_0_a, r_0_b, rdist_ij, rdist_jk, rij, rjk
4193 REAL(
dp),
DIMENSION(3) :: ftmp_i, ftmp_j, ftmp_k, ss, xij, xjk, &
4195 TYPE(particle_list_type),
POINTER :: particles_i
4196 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
4201 NULLIFY (particles_i)
4203 IF (
PRESENT(particles))
THEN
4204 my_particles => particles
4206 cpassert(
PRESENT(subsys))
4208 my_particles => particles_i%els
4210 n_atoms_to_a = colvar%coord_param%n_atoms_to
4211 n_atoms_to_b = colvar%coord_param%n_atoms_to_b
4212 n_atoms_from = colvar%coord_param%n_atoms_from
4213 p_a = colvar%coord_param%nncrd
4214 q_a = colvar%coord_param%ndcrd
4215 r_0_a = colvar%coord_param%r_0
4216 p_b = colvar%coord_param%nncrd_b
4217 q_b = colvar%coord_param%ndcrd_b
4218 r_0_b = colvar%coord_param%r_0_b
4221 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
4222 DO ii = 1, n_atoms_from
4223 i = colvar%coord_param%i_at_from(ii)
4224 CALL get_coordinates(colvar, i, xpi, my_particles)
4225 DO jj = 1, n_atoms_to_a
4226 j = colvar%coord_param%i_at_to(jj)
4227 CALL get_coordinates(colvar, j, xpj, my_particles)
4230 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
4232 xij = matmul(cell%hmat, ss)
4233 rij = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
4234 IF (rij < 1.0e-8_dp) cycle
4235 rdist_ij = rij/r_0_a
4236 IF (abs(1.0_dp - rdist_ij) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4237 num_ij = (1.0_dp - rdist_ij**p_a)
4238 invden_ij = 1.0_dp/(1.0_dp - rdist_ij**q_a)
4239 func_ij = num_ij*invden_ij
4240 IF (rij < 1.0e-8_dp)
THEN
4244 dfunc_ij = (-p_a*rdist_ij**(p_a - 1)*invden_ij &
4245 + num_ij*(invden_ij)**2*q_a*rdist_ij**(q_a - 1))/(rij*r_0_a)
4249 func_ij = real(p_a, kind=
dp)/real(q_a, kind=
dp)
4250 dfunc_ij = real(p_a, kind=
dp)*real((-q_a + p_a), kind=
dp)/(real(2*q_a, kind=
dp)*r_0_a)
4252 IF (n_atoms_to_b /= 0)
THEN
4254 DO kk = 1, n_atoms_to_b
4255 k = colvar%coord_param%i_at_to_b(kk)
4257 CALL get_coordinates(colvar, k, xpk, my_particles)
4258 ss = matmul(cell%h_inv, xpj(:) - xpk(:))
4260 xjk = matmul(cell%hmat, ss)
4261 rjk = sqrt(xjk(1)**2 + xjk(2)**2 + xjk(3)**2)
4262 IF (rjk < 1.0e-8_dp) cycle
4263 rdist_jk = rjk/r_0_b
4264 IF (abs(1.0_dp - rdist_jk) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4265 num_jk = (1.0_dp - rdist_jk**p_b)
4266 invden_jk = 1.0_dp/(1.0_dp - rdist_jk**q_b)
4267 func_jk = num_jk*invden_jk
4268 IF (rjk < 1.0e-8_dp)
THEN
4272 dfunc_jk = (-p_b*rdist_jk**(p_b - 1)*invden_jk &
4273 + num_jk*(invden_jk)**2*q_b*rdist_jk**(q_b - 1))/(rjk*r_0_b)
4277 func_jk = real(p_b, kind=
dp)/real(q_b, kind=
dp)
4278 dfunc_jk = real(p_b, kind=
dp)*real((-q_b + p_b), kind=
dp)/(real(2*q_b, kind=
dp)*r_0_b)
4280 func_k = func_k + func_jk
4281 ftmp_k = -func_ij*dfunc_jk*xjk
4282 CALL put_derivative(colvar, n_atoms_from + n_atoms_to_a + kk, ftmp_k)
4284 ftmp_j = -dfunc_ij*xij*func_jk + func_ij*dfunc_jk*xjk
4285 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4290 ftmp_j = -dfunc_ij*xij
4291 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4293 ncoord = ncoord + func_ij*func_k
4294 ftmp_i = dfunc_ij*xij*func_k
4295 CALL put_derivative(colvar, ii, ftmp_i)
4298 colvar%ss = ncoord*inv_n_atoms_from
4299 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
4300 END SUBROUTINE coord_colvar
4309 SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
4311 TYPE(colvar_type),
POINTER :: colvar
4312 TYPE(cell_type),
POINTER :: cell
4313 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
4314 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
4315 POINTER :: particles
4317 INTEGER :: i, ii, j, jj, n_coord_from, n_coord_to, &
4319 REAL(
dp) :: den_n, den_q, fscalar, ftemp_i(3), inv_den_n, inv_den_q, lambda, num_n, num_q, &
4320 qfunc, r12, r_cut, rfact, rij(3), rpi(3), rpj(3)
4321 REAL(
dp),
DIMENSION(:),
POINTER :: dqfunc_dnl, expnl, nlcoord, sum_rij
4322 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: dnlcoord, dqfunc_dr
4323 TYPE(particle_list_type),
POINTER :: particles_i
4324 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
4329 NULLIFY (particles_i)
4331 IF (
PRESENT(particles))
THEN
4332 my_particles => particles
4334 cpassert(
PRESENT(subsys))
4336 my_particles => particles_i%els
4339 n_dist_from = colvar%mindist_param%n_dist_from
4340 n_coord_from = colvar%mindist_param%n_coord_from
4341 n_coord_to = colvar%mindist_param%n_coord_to
4342 p = colvar%mindist_param%p_exp
4343 q = colvar%mindist_param%q_exp
4344 r_cut = colvar%mindist_param%r_cut
4345 lambda = colvar%mindist_param%lambda
4347 NULLIFY (nlcoord, dnlcoord, dqfunc_dr, dqfunc_dnl, expnl, sum_rij)
4348 ALLOCATE (nlcoord(n_coord_from))
4349 ALLOCATE (dnlcoord(3, n_coord_from, n_coord_to))
4350 ALLOCATE (expnl(n_coord_from))
4351 ALLOCATE (sum_rij(n_coord_from))
4352 ALLOCATE (dqfunc_dr(3, n_dist_from, n_coord_from))
4353 ALLOCATE (dqfunc_dnl(n_coord_from))
4360 DO i = 1, n_coord_from
4361 ii = colvar%mindist_param%i_coord_from(i)
4362 rpi = my_particles(ii)%r(1:3)
4363 DO j = 1, n_coord_to
4364 jj = colvar%mindist_param%i_coord_to(j)
4365 rpj = my_particles(jj)%r(1:3)
4366 rij =
pbc(rpj, rpi, cell)
4367 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4369 num_n = 1.0_dp - rfact**p
4370 den_n = 1.0_dp - rfact**q
4371 inv_den_n = 1.0_dp/den_n
4372 IF (abs(inv_den_n) < 1.e-10_dp)
THEN
4373 inv_den_n = 1.e-10_dp
4377 fscalar = (-p*rfact**(p - 1) + num_n*q*rfact**(q - 1)*inv_den_n)*inv_den_n/(r_cut*r12)
4379 dnlcoord(1, i, j) = rij(1)*fscalar
4380 dnlcoord(2, i, j) = rij(2)*fscalar
4381 dnlcoord(3, i, j) = rij(3)*fscalar
4383 nlcoord(i) = nlcoord(i) + num_n*inv_den_n
4385 expnl(i) = exp(lambda*nlcoord(i))
4389 den_q = den_q + expnl(i)
4391 inv_den_q = 1.0_dp/den_q
4398 DO i = 1, n_dist_from
4399 ii = colvar%mindist_param%i_dist_from(i)
4400 rpi = my_particles(ii)%r(1:3)
4401 DO j = 1, n_coord_from
4402 jj = colvar%mindist_param%i_coord_from(j)
4403 rpj = my_particles(jj)%r(1:3)
4404 rij =
pbc(rpj, rpi, cell)
4405 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4410 num_q = num_q + r12*expnl(j)
4412 sum_rij(j) = sum_rij(j) + r12
4413 dqfunc_dr(1, i, j) = expnl(j)*rij(1)/r12
4414 dqfunc_dr(2, i, j) = expnl(j)*rij(2)/r12
4415 dqfunc_dr(3, i, j) = expnl(j)*rij(3)/r12
4422 qfunc = num_q*inv_den_q
4423 dqfunc_dr = dqfunc_dr*inv_den_q
4430 DO i = 1, n_coord_from
4431 dqfunc_dnl(i) = lambda*expnl(i)*inv_den_q*(sum_rij(i) - num_q*inv_den_q)
4435 DO i = 1, n_dist_from
4436 DO j = 1, n_coord_from
4437 ftemp_i(1) = dqfunc_dr(1, i, j)
4438 ftemp_i(2) = dqfunc_dr(2, i, j)
4439 ftemp_i(3) = dqfunc_dr(3, i, j)
4441 CALL put_derivative(colvar, i, ftemp_i)
4442 CALL put_derivative(colvar, j + n_dist_from, -ftemp_i)
4446 DO i = 1, n_coord_from
4447 DO j = 1, n_coord_to
4448 ftemp_i(1) = dqfunc_dnl(i)*dnlcoord(1, i, j)
4449 ftemp_i(2) = dqfunc_dnl(i)*dnlcoord(2, i, j)
4450 ftemp_i(3) = dqfunc_dnl(i)*dnlcoord(3, i, j)
4452 CALL put_derivative(colvar, i + n_dist_from, ftemp_i)
4453 CALL put_derivative(colvar, j + n_dist_from + n_coord_from, -ftemp_i)
4458 DEALLOCATE (nlcoord)
4459 DEALLOCATE (dnlcoord)
4461 DEALLOCATE (dqfunc_dr)
4462 DEALLOCATE (sum_rij)
4463 DEALLOCATE (dqfunc_dnl)
4465 END SUBROUTINE mindist_colvar
4475 SUBROUTINE combine_colvar(colvar, cell, subsys, particles)
4476 TYPE(colvar_type),
POINTER :: colvar
4477 TYPE(cell_type),
POINTER :: cell
4478 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
4479 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
4480 POINTER :: particles
4482 CHARACTER(LEN=default_string_length) :: def_error, this_error
4483 CHARACTER(LEN=default_string_length), &
4484 ALLOCATABLE,
DIMENSION(:) :: my_par
4485 INTEGER :: i, ii, j, ncolv, ndim
4487 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dss_vals, my_val, ss_vals
4488 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi
4489 TYPE(particle_list_type),
POINTER :: particles_i
4490 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
4493 IF (
PRESENT(particles))
THEN
4494 my_particles => particles
4496 cpassert(
PRESENT(subsys))
4498 my_particles => particles_i%els
4501 ncolv =
SIZE(colvar%combine_cvs_param%colvar_p)
4502 ALLOCATE (ss_vals(ncolv))
4503 ALLOCATE (dss_vals(ncolv))
4507 CALL colvar_recursive_eval(colvar%combine_cvs_param%colvar_p(i)%colvar, cell, my_particles)
4508 ss_vals(i) = colvar%combine_cvs_param%colvar_p(i)%colvar%ss
4513 ndim =
SIZE(colvar%combine_cvs_param%c_parameters) + &
4514 SIZE(colvar%combine_cvs_param%variables)
4515 ALLOCATE (my_par(ndim))
4516 my_par(1:
SIZE(colvar%combine_cvs_param%variables)) = colvar%combine_cvs_param%variables
4517 my_par(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%c_parameters
4518 ALLOCATE (my_val(ndim))
4519 my_val(1:
SIZE(colvar%combine_cvs_param%variables)) = ss_vals
4520 my_val(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%v_parameters
4521 CALL parsef(1, trim(colvar%combine_cvs_param%function), my_par)
4522 colvar%ss =
evalf(1, my_val)
4524 dss_vals(i) =
evalfd(1, i, my_val, colvar%combine_cvs_param%dx, err)
4525 IF ((abs(err) > colvar%combine_cvs_param%lerr))
THEN
4526 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
4527 WRITE (def_error,
"(A,G12.6,A)")
"(", colvar%combine_cvs_param%lerr,
")"
4530 CALL cp_warn(__location__, &
4531 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
4532 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
4533 trim(def_error)//
' . ')
4541 ALLOCATE (fi(3, colvar%n_atom_s))
4544 DO j = 1, colvar%combine_cvs_param%colvar_p(i)%colvar%n_atom_s
4546 fi(:, ii) = colvar%combine_cvs_param%colvar_p(i)%colvar%dsdr(:, j)*dss_vals(i)
4550 DO i = 1, colvar%n_atom_s
4551 CALL put_derivative(colvar, i, fi(:, i))
4555 DEALLOCATE (ss_vals)
4556 DEALLOCATE (dss_vals)
4557 END SUBROUTINE combine_colvar
4573 SUBROUTINE reaction_path_colvar(colvar, cell, subsys, particles)
4574 TYPE(colvar_type),
POINTER :: colvar
4575 TYPE(cell_type),
POINTER :: cell
4576 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
4577 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
4578 POINTER :: particles
4580 TYPE(particle_list_type),
POINTER :: particles_i
4581 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
4584 IF (
PRESENT(particles))
THEN
4585 my_particles => particles
4587 cpassert(
PRESENT(subsys))
4589 my_particles => particles_i%els
4592 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4593 CALL rpath_dist_rmsd(colvar, my_particles)
4594 ELSEIF (colvar%reaction_path_param%rmsd)
THEN
4595 CALL rpath_rmsd(colvar, my_particles)
4597 CALL rpath_colvar(colvar, cell, my_particles)
4600 END SUBROUTINE reaction_path_colvar
4611 SUBROUTINE rpath_colvar(colvar, cell, particles)
4612 TYPE(colvar_type),
POINTER :: colvar
4613 TYPE(cell_type),
POINTER :: cell
4614 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
4616 INTEGER :: i, iend, ii, istart, j, k, ncolv, nconf
4617 REAL(
dp) :: lambda, step_size
4618 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: s1, ss_vals
4619 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, f_vals, fi, s1v
4620 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
4622 istart = colvar%reaction_path_param%function_bounds(1)
4623 iend = colvar%reaction_path_param%function_bounds(2)
4625 nconf = colvar%reaction_path_param%nr_frames
4626 step_size = colvar%reaction_path_param%step_size
4627 ncolv = colvar%reaction_path_param%n_components
4628 lambda = colvar%reaction_path_param%lambda
4629 ALLOCATE (f_vals(ncolv, istart:iend))
4630 f_vals(:, :) = colvar%reaction_path_param%f_vals
4631 ALLOCATE (ss_vals(ncolv))
4634 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4635 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4638 ALLOCATE (s1v(2, istart:iend))
4639 ALLOCATE (ds1v(ncolv, 2, istart:iend))
4642 ALLOCATE (ds1(ncolv, 2))
4645 s1v(1, k) = real(k, kind=
dp)*step_size*exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4646 s1v(2, k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4648 ds1v(j, 1, k) = f_vals(j, k)*s1v(1, k)
4649 ds1v(j, 2, k) = f_vals(j, k)*s1v(2, k)
4653 s1(i) = accurate_sum(s1v(i, :))
4655 ds1(j, i) = accurate_sum(ds1v(j, i, :))
4659 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4661 ALLOCATE (fi(3, colvar%n_atom_s))
4665 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
4667 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)*lambda* &
4668 (ds1(i, 1)/s1(2)/real(nconf - 1,
dp) - colvar%ss*ds1(i, 2)/s1(2))*2.0_dp
4672 DO i = 1, colvar%n_atom_s
4673 CALL put_derivative(colvar, i, fi(:, i))
4678 DEALLOCATE (ss_vals)
4684 END SUBROUTINE rpath_colvar
4695 SUBROUTINE rpath_dist_rmsd(colvar, particles)
4696 TYPE(colvar_type),
POINTER :: colvar
4697 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
4699 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4700 INTEGER,
DIMENSION(:),
POINTER :: iatom
4701 REAL(
dp) :: lambda, my_rmsd, s1(2), sum_exp
4702 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, vec_dif
4703 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dvec_dif, fi, riat, s1v
4704 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4705 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4706 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4708 nconf = colvar%reaction_path_param%nr_frames
4709 rmsd_atom = colvar%reaction_path_param%n_components
4710 lambda = colvar%reaction_path_param%lambda
4711 path_conf => colvar%reaction_path_param%r_ref
4712 iatom => colvar%reaction_path_param%i_rmsd
4714 natom =
SIZE(particles)
4716 ALLOCATE (r0(3*natom))
4717 ALLOCATE (r(3*natom))
4718 ALLOCATE (riat(3, rmsd_atom))
4719 ALLOCATE (vec_dif(rmsd_atom))
4720 ALLOCATE (dvec_dif(3, rmsd_atom))
4721 ALLOCATE (s1v(2, nconf))
4722 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4723 ALLOCATE (ds1(3, rmsd_atom, 2))
4726 r0(ii + 1) = particles(i)%r(1)
4727 r0(ii + 2) = particles(i)%r(2)
4728 r0(ii + 3) = particles(i)%r(3)
4731 DO iat = 1, rmsd_atom
4733 riat(:, iat) = particles(ii)%r
4739 r(ii + 1) = path_conf(ii + 1, ik)
4740 r(ii + 2) = path_conf(ii + 2, ik)
4741 r(ii + 3) = path_conf(ii + 3, ik)
4744 CALL rmsd3(particles, r, r0, output_unit=-1, my_val=my_rmsd, rotate=.true.)
4747 DO iat = 1, rmsd_atom
4750 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 &
4751 + (riat(3, iat) - r(ii + 3))**2
4752 sum_exp = sum_exp + vec_dif(iat)
4755 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*sum_exp)
4756 s1v(2, ik) = exp(-lambda*sum_exp)
4757 DO iat = 1, rmsd_atom
4760 ds1v(1, iat, 1, ik) = r(ii + 1)*s1v(1, ik)
4761 ds1v(1, iat, 2, ik) = r(ii + 1)*s1v(2, ik)
4762 ds1v(2, iat, 1, ik) = r(ii + 2)*s1v(1, ik)
4763 ds1v(2, iat, 2, ik) = r(ii + 2)*s1v(2, ik)
4764 ds1v(3, iat, 1, ik) = r(ii + 3)*s1v(1, ik)
4765 ds1v(3, iat, 2, ik) = r(ii + 3)*s1v(2, ik)
4769 s1(1) = accurate_sum(s1v(1, :))
4770 s1(2) = accurate_sum(s1v(2, :))
4772 DO iat = 1, rmsd_atom
4773 ds1(1, iat, i) = accurate_sum(ds1v(1, iat, i, :))
4774 ds1(2, iat, i) = accurate_sum(ds1v(2, iat, i, :))
4775 ds1(3, iat, i) = accurate_sum(ds1v(3, iat, i, :))
4779 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4781 ALLOCATE (fi(3, rmsd_atom))
4783 DO iat = 1, rmsd_atom
4784 fi(1, iat) = 2.0_dp*lambda/s1(2)/real(nconf - 1,
dp)*(ds1(1, iat, 1) - ds1(1, iat, 2)*s1(1)/s1(2))
4785 fi(2, iat) = 2.0_dp*lambda/s1(2)/real(nconf - 1,
dp)*(ds1(2, iat, 1) - ds1(2, iat, 2)*s1(1)/s1(2))
4786 fi(3, iat) = 2.0_dp*lambda/s1(2)/real(nconf - 1,
dp)*(ds1(3, iat, 1) - ds1(3, iat, 2)*s1(1)/s1(2))
4787 CALL put_derivative(colvar, iat, fi(:, iat))
4794 DEALLOCATE (vec_dif)
4795 DEALLOCATE (dvec_dif)
4800 END SUBROUTINE rpath_dist_rmsd
4807 SUBROUTINE rpath_rmsd(colvar, particles)
4808 TYPE(colvar_type),
POINTER :: colvar
4809 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
4811 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4812 INTEGER,
DIMENSION(:),
POINTER :: iatom
4813 REAL(
dp) :: lambda, my_rmsd, s1(2)
4814 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0
4815 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi, riat, s1v
4816 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4817 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4818 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4819 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
4820 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
4822 nconf = colvar%reaction_path_param%nr_frames
4823 rmsd_atom = colvar%reaction_path_param%n_components
4824 lambda = colvar%reaction_path_param%lambda
4825 path_conf => colvar%reaction_path_param%r_ref
4826 iatom => colvar%reaction_path_param%i_rmsd
4828 natom =
SIZE(particles)
4830 ALLOCATE (r0(3*natom))
4831 ALLOCATE (r(3*natom))
4832 ALLOCATE (riat(3, rmsd_atom))
4833 ALLOCATE (s1v(2, nconf))
4834 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4835 ALLOCATE (ds1(3, rmsd_atom, 2))
4836 ALLOCATE (drmsd(3, natom))
4838 ALLOCATE (weight(natom))
4842 r0(ii + 1) = particles(i)%r(1)
4843 r0(ii + 2) = particles(i)%r(2)
4844 r0(ii + 3) = particles(i)%r(3)
4847 DO iat = 1, rmsd_atom
4849 riat(:, iat) = particles(ii)%r
4854 DO iat = 1, rmsd_atom
4862 r(ii + 1) = path_conf(ii + 1, ik)
4863 r(ii + 2) = path_conf(ii + 2, ik)
4864 r(ii + 3) = path_conf(ii + 3, ik)
4867 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
4868 rotate=.false., drmsd3=drmsd)
4870 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*my_rmsd)
4871 s1v(2, ik) = exp(-lambda*my_rmsd)
4872 DO iat = 1, rmsd_atom
4874 ds1v(1, iat, 1, ik) = drmsd(1, i)*s1v(1, ik)
4875 ds1v(1, iat, 2, ik) = drmsd(1, i)*s1v(2, ik)
4876 ds1v(2, iat, 1, ik) = drmsd(2, i)*s1v(1, ik)
4877 ds1v(2, iat, 2, ik) = drmsd(2, i)*s1v(2, ik)
4878 ds1v(3, iat, 1, ik) = drmsd(3, i)*s1v(1, ik)
4879 ds1v(3, iat, 2, ik) = drmsd(3, i)*s1v(2, ik)
4883 s1(1) = accurate_sum(s1v(1, :))
4884 s1(2) = accurate_sum(s1v(2, :))
4886 DO iat = 1, rmsd_atom
4887 ds1(1, iat, i) = accurate_sum(ds1v(1, iat, i, :))
4888 ds1(2, iat, i) = accurate_sum(ds1v(2, iat, i, :))
4889 ds1(3, iat, i) = accurate_sum(ds1v(3, iat, i, :))
4893 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4895 ALLOCATE (fi(3, rmsd_atom))
4897 DO iat = 1, rmsd_atom
4898 fi(1, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(1, iat, 1) - ds1(1, iat, 2)*s1(1)/s1(2))
4899 fi(2, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(2, iat, 1) - ds1(2, iat, 2)*s1(1)/s1(2))
4900 fi(3, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(3, iat, 1) - ds1(3, iat, 2)*s1(1)/s1(2))
4901 CALL put_derivative(colvar, iat, fi(:, iat))
4914 END SUBROUTINE rpath_rmsd
4926 SUBROUTINE distance_from_path_colvar(colvar, cell, subsys, particles)
4927 TYPE(colvar_type),
POINTER :: colvar
4928 TYPE(cell_type),
POINTER :: cell
4929 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
4930 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
4931 POINTER :: particles
4933 TYPE(particle_list_type),
POINTER :: particles_i
4934 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
4937 IF (
PRESENT(particles))
THEN
4938 my_particles => particles
4940 cpassert(
PRESENT(subsys))
4942 my_particles => particles_i%els
4945 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4946 CALL dpath_dist_rmsd(colvar, my_particles)
4947 ELSEIF (colvar%reaction_path_param%rmsd)
THEN
4948 CALL dpath_rmsd(colvar, my_particles)
4950 CALL dpath_colvar(colvar, cell, my_particles)
4953 END SUBROUTINE distance_from_path_colvar
4965 SUBROUTINE dpath_colvar(colvar, cell, particles)
4966 TYPE(colvar_type),
POINTER :: colvar
4967 TYPE(cell_type),
POINTER :: cell
4968 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
4970 INTEGER :: i, iend, ii, istart, j, k, ncolv
4971 REAL(
dp) :: lambda, s1
4972 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: ds1, s1v, ss_vals
4973 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1v, f_vals, fi
4975 istart = colvar%reaction_path_param%function_bounds(1)
4976 iend = colvar%reaction_path_param%function_bounds(2)
4978 ncolv = colvar%reaction_path_param%n_components
4979 lambda = colvar%reaction_path_param%lambda
4980 ALLOCATE (f_vals(ncolv, istart:iend))
4981 f_vals(:, :) = colvar%reaction_path_param%f_vals
4982 ALLOCATE (ss_vals(ncolv))
4985 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4986 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4989 ALLOCATE (s1v(istart:iend))
4990 ALLOCATE (ds1v(ncolv, istart:iend))
4991 ALLOCATE (ds1(ncolv))
4994 s1v(k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4996 ds1v(j, k) = f_vals(j, k)*s1v(k)
5000 s1 = accurate_sum(s1v(:))
5002 ds1(j) = accurate_sum(ds1v(j, :))
5004 colvar%ss = -1.0_dp/lambda*log(s1)
5006 ALLOCATE (fi(3, colvar%n_atom_s))
5010 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
5012 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)* &
5013 2.0_dp*(ss_vals(i) - ds1(i)/s1)
5017 DO i = 1, colvar%n_atom_s
5018 CALL put_derivative(colvar, i, fi(:, i))
5023 DEALLOCATE (ss_vals)
5028 END SUBROUTINE dpath_colvar
5039 SUBROUTINE dpath_dist_rmsd(colvar, particles)
5041 TYPE(colvar_type),
POINTER :: colvar
5042 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
5044 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5045 INTEGER,
DIMENSION(:),
POINTER :: iatom
5046 REAL(
dp) :: lambda, s1, sum_exp
5047 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v, vec_dif
5048 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, dvec_dif, fi, riat
5049 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5050 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5052 nconf = colvar%reaction_path_param%nr_frames
5053 rmsd_atom = colvar%reaction_path_param%n_components
5054 lambda = colvar%reaction_path_param%lambda
5055 path_conf => colvar%reaction_path_param%r_ref
5056 iatom => colvar%reaction_path_param%i_rmsd
5058 natom =
SIZE(particles)
5060 ALLOCATE (r0(3*natom))
5061 ALLOCATE (r(3*natom))
5062 ALLOCATE (riat(3, rmsd_atom))
5063 ALLOCATE (vec_dif(rmsd_atom))
5064 ALLOCATE (dvec_dif(3, rmsd_atom))
5065 ALLOCATE (s1v(nconf))
5066 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5067 ALLOCATE (ds1(3, rmsd_atom))
5070 r0(ii + 1) = particles(i)%r(1)
5071 r0(ii + 2) = particles(i)%r(2)
5072 r0(ii + 3) = particles(i)%r(3)
5075 DO iat = 1, rmsd_atom
5077 riat(:, iat) = particles(ii)%r
5083 r(ii + 1) = path_conf(ii + 1, ik)
5084 r(ii + 2) = path_conf(ii + 2, ik)
5085 r(ii + 3) = path_conf(ii + 3, ik)
5088 CALL rmsd3(particles, r, r0, output_unit=-1, rotate=.true.)
5091 DO iat = 1, rmsd_atom
5094 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 + (riat(3, iat) - r(ii + 3))**2
5095 sum_exp = sum_exp + vec_dif(iat)
5096 dvec_dif(1, iat) = r(ii + 1)
5097 dvec_dif(2, iat) = r(ii + 2)
5098 dvec_dif(3, iat) = r(ii + 3)
5100 s1v(ik) = exp(-lambda*sum_exp)
5101 DO iat = 1, rmsd_atom
5102 ds1v(1, iat, ik) = dvec_dif(1, iat)*s1v(ik)
5103 ds1v(2, iat, ik) = dvec_dif(2, iat)*s1v(ik)
5104 ds1v(3, iat, ik) = dvec_dif(3, iat)*s1v(ik)
5108 s1 = accurate_sum(s1v(:))
5109 DO iat = 1, rmsd_atom
5110 ds1(1, iat) = accurate_sum(ds1v(1, iat, :))
5111 ds1(2, iat) = accurate_sum(ds1v(2, iat, :))
5112 ds1(3, iat) = accurate_sum(ds1v(3, iat, :))
5114 colvar%ss = -1.0_dp/lambda*log(s1)
5116 ALLOCATE (fi(3, rmsd_atom))
5118 DO iat = 1, rmsd_atom
5119 fi(:, iat) = 2.0_dp*(riat(:, iat) - ds1(:, iat)/s1)
5120 CALL put_derivative(colvar, iat, fi(:, iat))
5127 DEALLOCATE (vec_dif)
5128 DEALLOCATE (dvec_dif)
5132 END SUBROUTINE dpath_dist_rmsd
5139 SUBROUTINE dpath_rmsd(colvar, particles)
5141 TYPE(colvar_type),
POINTER :: colvar
5142 TYPE(particle_type),
DIMENSION(:),
POINTER :: particles
5144 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5145 INTEGER,
DIMENSION(:),
POINTER :: iatom
5146 REAL(
dp) :: lambda, my_rmsd, s1
5147 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v
5148 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, fi, riat
5149 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5150 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5151 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
5152 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
5154 nconf = colvar%reaction_path_param%nr_frames
5155 rmsd_atom = colvar%reaction_path_param%n_components
5156 lambda = colvar%reaction_path_param%lambda
5157 path_conf => colvar%reaction_path_param%r_ref
5158 iatom => colvar%reaction_path_param%i_rmsd
5160 natom =
SIZE(particles)
5162 ALLOCATE (r0(3*natom))
5163 ALLOCATE (r(3*natom))
5164 ALLOCATE (riat(3, rmsd_atom))
5165 ALLOCATE (s1v(nconf))
5166 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5167 ALLOCATE (ds1(3, rmsd_atom))
5168 ALLOCATE (drmsd(3, natom))
5170 ALLOCATE (weight(natom))
5174 r0(ii + 1) = particles(i)%r(1)
5175 r0(ii + 2) = particles(i)%r(2)
5176 r0(ii + 3) = particles(i)%r(3)
5179 DO iat = 1, rmsd_atom
5181 riat(:, iat) = particles(ii)%r
5186 DO iat = 1, rmsd_atom
5194 r(ii + 1) = path_conf(ii + 1, ik)
5195 r(ii + 2) = path_conf(ii + 2, ik)
5196 r(ii + 3) = path_conf(ii + 3, ik)
5199 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
5200 rotate=.false., drmsd3=drmsd)
5202 s1v(ik) = exp(-lambda*my_rmsd)
5203 DO iat = 1, rmsd_atom
5205 ds1v(1, iat, ik) = drmsd(1, i)*s1v(ik)
5206 ds1v(2, iat, ik) = drmsd(2, i)*s1v(ik)
5207 ds1v(3, iat, ik) = drmsd(3, i)*s1v(ik)
5211 s1 = accurate_sum(s1v(:))
5212 DO iat = 1, rmsd_atom
5213 ds1(1, iat) = accurate_sum(ds1v(1, iat, :))
5214 ds1(2, iat) = accurate_sum(ds1v(2, iat, :))
5215 ds1(3, iat) = accurate_sum(ds1v(3, iat, :))
5217 colvar%ss = -1.0_dp/lambda*log(s1)
5219 ALLOCATE (fi(3, rmsd_atom))
5221 DO iat = 1, rmsd_atom
5222 fi(:, iat) = ds1(:, iat)/s1
5223 CALL put_derivative(colvar, iat, fi(:, iat))
5236 END SUBROUTINE dpath_rmsd
5247 SUBROUTINE population_colvar(colvar, cell, subsys, particles)
5248 TYPE(colvar_type),
POINTER :: colvar
5249 TYPE(cell_type),
POINTER :: cell
5250 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5251 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
5252 POINTER :: particles
5254 INTEGER :: i, ii, jj, n_atoms_from, n_atoms_to, &
5256 REAL(
dp) :: dfunc, dfunc_coord, ftmp(3), func, func_coord, inv_n_atoms_from, invden, n_0, &
5257 ncoord, norm, num, population, r12, r_0, rdist, sigma, ss(3), xij(3)
5258 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftmp_coord
5259 REAL(
dp),
DIMENSION(3) :: xpi, xpj
5260 TYPE(particle_list_type),
POINTER :: particles_i
5261 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
5266 NULLIFY (particles_i)
5268 IF (
PRESENT(particles))
THEN
5269 my_particles => particles
5271 cpassert(
PRESENT(subsys))
5273 my_particles => particles_i%els
5275 n_atoms_to = colvar%population_param%n_atoms_to
5276 n_atoms_from = colvar%population_param%n_atoms_from
5277 nncrd = colvar%population_param%nncrd
5278 ndcrd = colvar%population_param%ndcrd
5279 r_0 = colvar%population_param%r_0
5280 n_0 = colvar%population_param%n0
5281 sigma = colvar%population_param%sigma
5283 ALLOCATE (ftmp_coord(3, n_atoms_to))
5289 colvar%dsdr = 0.0_dp
5290 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
5292 norm = sqrt(
pi*2.0_dp)*sigma
5295 DO ii = 1, n_atoms_from
5296 i = colvar%population_param%i_at_from(ii)
5297 CALL get_coordinates(colvar, i, xpi, my_particles)
5298 DO jj = 1, n_atoms_to
5299 i = colvar%population_param%i_at_to(jj)
5300 CALL get_coordinates(colvar, i, xpj, my_particles)
5301 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
5303 xij = matmul(cell%hmat, ss)
5304 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
5305 IF (r12 < 1.0e-8_dp) cycle
5307 num = (1.0_dp - rdist**nncrd)
5308 invden = 1.0_dp/(1.0_dp - rdist**ndcrd)
5309 func_coord = num*invden
5310 dfunc_coord = (-nncrd*rdist**(nncrd - 1)*invden &
5311 + num*(invden)**2*ndcrd*rdist**(ndcrd - 1))/(r12*r_0)
5313 ncoord = ncoord + func_coord
5314 ftmp_coord(1, jj) = dfunc_coord*xij(1)
5315 ftmp_coord(2, jj) = dfunc_coord*xij(2)
5316 ftmp_coord(3, jj) = dfunc_coord*xij(3)
5319 func = exp(-(ncoord - n_0)**2/(2.0_dp*sigma*sigma))
5320 dfunc = -func*(ncoord - n_0)/(sigma*sigma)
5322 population = population + norm*func
5323 DO jj = 1, n_atoms_to
5324 ftmp(1) = ftmp_coord(1, jj)*dfunc
5325 ftmp(2) = ftmp_coord(2, jj)*dfunc
5326 ftmp(3) = ftmp_coord(3, jj)*dfunc
5327 CALL put_derivative(colvar, ii, ftmp)
5328 ftmp(1) = -ftmp_coord(1, jj)*dfunc
5329 ftmp(2) = -ftmp_coord(2, jj)*dfunc
5330 ftmp(3) = -ftmp_coord(3, jj)*dfunc
5331 CALL put_derivative(colvar, n_atoms_from + jj, ftmp)
5335 colvar%ss = population
5336 END SUBROUTINE population_colvar
5348 SUBROUTINE gyration_radius_colvar(colvar, cell, subsys, particles)
5350 TYPE(colvar_type),
POINTER :: colvar
5351 TYPE(cell_type),
POINTER :: cell
5352 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5353 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
5354 POINTER :: particles
5356 INTEGER :: i, ii, n_atoms
5357 REAL(
dp) :: dri2, func, gyration, inv_n, mass_tot, mi
5358 REAL(
dp),
DIMENSION(3) :: dfunc, dxi, ftmp, ss, xpcom, xpi
5359 TYPE(particle_list_type),
POINTER :: particles_i
5360 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
5362 NULLIFY (particles_i, my_particles)
5364 IF (
PRESENT(particles))
THEN
5365 my_particles => particles
5367 cpassert(
PRESENT(subsys))
5369 my_particles => particles_i%els
5371 n_atoms = colvar%gyration_param%n_atoms
5372 inv_n = 1.0_dp/n_atoms
5378 i = colvar%gyration_param%i_at(ii)
5379 CALL get_coordinates(colvar, i, xpi, my_particles)
5380 CALL get_mass(colvar, i, mi, my_particles)
5381 xpcom(:) = xpcom(:) + xpi(:)*mi
5382 mass_tot = mass_tot + mi
5384 xpcom(:) = xpcom(:)/mass_tot
5390 i = colvar%gyration_param%i_at(ii)
5391 CALL get_coordinates(colvar, i, xpi, my_particles)
5392 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5394 dxi = matmul(cell%hmat, ss)
5395 dri2 = (dxi(1)**2 + dxi(2)**2 + dxi(3)**2)
5397 dfunc(:) = dfunc(:) + dxi(:)
5399 gyration = sqrt(inv_n*func)
5402 i = colvar%gyration_param%i_at(ii)
5403 CALL get_coordinates(colvar, i, xpi, my_particles)
5404 CALL get_mass(colvar, i, mi, my_particles)
5405 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5407 dxi = matmul(cell%hmat, ss)
5408 ftmp(1) = dxi(1) - dfunc(1)*mi/mass_tot
5409 ftmp(2) = dxi(2) - dfunc(2)*mi/mass_tot
5410 ftmp(3) = dxi(3) - dfunc(3)*mi/mass_tot
5411 ftmp(:) = ftmp(:)*inv_n/gyration
5412 CALL put_derivative(colvar, ii, ftmp)
5414 colvar%ss = gyration
5416 END SUBROUTINE gyration_radius_colvar
5427 SUBROUTINE rmsd_colvar(colvar, subsys, particles)
5428 TYPE(colvar_type),
POINTER :: colvar
5429 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5430 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
5431 POINTER :: particles
5433 CALL rmsd_colvar_low(colvar, subsys, particles)
5434 END SUBROUTINE rmsd_colvar
5449 SUBROUTINE rmsd_colvar_low(colvar, subsys, particles)
5451 TYPE(colvar_type),
POINTER :: colvar
5452 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5453 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
5454 POINTER :: particles
5456 INTEGER :: i, ii, natom, nframes
5457 REAL(kind=
dp) :: cv_val, f1, ftmp(3)
5458 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: der, r,
rmsd
5459 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r0
5460 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: drmsd
5461 REAL(kind=
dp),
DIMENSION(:),
POINTER :: weights
5462 TYPE(particle_list_type),
POINTER :: particles_i
5463 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
5465 NULLIFY (my_particles, particles_i, weights)
5467 IF (
PRESENT(particles))
THEN
5468 my_particles => particles
5470 cpassert(
PRESENT(subsys))
5472 my_particles => particles_i%els
5475 natom =
SIZE(my_particles)
5476 nframes = colvar%rmsd_param%nr_frames
5477 ALLOCATE (drmsd(3, natom, nframes))
5480 ALLOCATE (r0(3*natom, nframes))
5481 ALLOCATE (
rmsd(nframes))
5482 ALLOCATE (der(nframes))
5483 ALLOCATE (r(3*natom))
5485 weights => colvar%rmsd_param%weights
5488 r(ii + 1) = my_particles(i)%r(1)
5489 r(ii + 2) = my_particles(i)%r(2)
5490 r(ii + 3) = my_particles(i)%r(3)
5492 r0(:, :) = colvar%rmsd_param%r_ref
5495 CALL rmsd3(my_particles, r, r0(:, 1), output_unit=-1, weights=weights, my_val=
rmsd(1), rotate=.false., drmsd3=drmsd(:, :, 1))
5497 IF (nframes == 2)
THEN
5498 CALL rmsd3(my_particles, r, r0(:, 2), output_unit=-1, weights=weights, &
5499 my_val=
rmsd(2), rotate=.false., drmsd3=drmsd(:, :, 2))
5505 der(1) = f1 - cv_val*f1
5507 der(2) = -f1 - cv_val*f1
5509 DO i = 1, colvar%rmsd_param%n_atoms
5510 ii = colvar%rmsd_param%i_rmsd(i)
5511 IF (weights(ii) > 0.0_dp)
THEN
5512 ftmp(1) = der(1)*drmsd(1, ii, 1) + der(2)*drmsd(1, ii, 2)
5513 ftmp(2) = der(1)*drmsd(2, ii, 1) + der(2)*drmsd(2, ii, 2)
5514 ftmp(3) = der(1)*drmsd(3, ii, 1) + der(2)*drmsd(3, ii, 2)
5515 CALL put_derivative(colvar, i, ftmp)
5518 ELSE IF (nframes == 1)
THEN
5521 cv_val = sqrt(
rmsd(1))
5523 IF (cv_val /= 0.0_dp) f1 = 0.5_dp/cv_val
5524 DO i = 1, colvar%rmsd_param%n_atoms
5525 ii = colvar%rmsd_param%i_rmsd(i)
5526 IF (weights(ii) > 0.0_dp)
THEN
5527 ftmp(1) = f1*drmsd(1, ii, 1)
5528 ftmp(2) = f1*drmsd(2, ii, 1)
5529 ftmp(3) = f1*drmsd(3, ii, 1)
5530 CALL put_derivative(colvar, i, ftmp)
5534 cpabort(
"RMSD implemented only for 1 and 2 reference frames!")
5544 END SUBROUTINE rmsd_colvar_low
5556 SUBROUTINE ring_puckering_colvar(colvar, cell, subsys, particles)
5557 TYPE(colvar_type),
POINTER :: colvar
5558 TYPE(cell_type),
POINTER :: cell
5559 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5560 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
5561 POINTER :: particles
5563 INTEGER :: i, ii, j, jj, m, nring
5564 REAL(kind=
dp) :: a, at, b, da, db, ds, kr, rpxpp, svar
5565 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cosj, sinj, z
5566 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r
5567 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: nforce, zforce
5568 REAL(kind=
dp),
DIMENSION(3) :: ftmp, nv, r0, rp, rpp, uv
5569 REAL(kind=
dp),
DIMENSION(3, 3) :: dnvp, dnvpp
5570 TYPE(particle_list_type),
POINTER :: particles_i
5571 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
5574 IF (
PRESENT(particles))
THEN
5575 my_particles => particles
5577 cpassert(
PRESENT(subsys))
5579 my_particles => particles_i%els
5582 nring = colvar%ring_puckering_param%nring
5583 ALLOCATE (r(3, nring), z(nring), cosj(nring), sinj(nring))
5584 ALLOCATE (nforce(3, 3, nring), zforce(nring, nring, 3))
5586 i = colvar%ring_puckering_param%atoms(ii)
5587 CALL get_coordinates(colvar, i, r(:, ii), my_particles)
5592 r(:, ii) =
pbc(r(:, ii), r0, cell)
5597 r0(:) = r0(:) + r(:, ii)
5599 kr = 1._dp/real(nring, kind=
dp)
5602 r(:, ii) = r(:, ii) - r0(:)
5608 cosj(ii) = cos(
twopi*(ii - 1)*kr)
5609 sinj(ii) = sin(
twopi*(ii - 1)*kr)
5610 rp(:) = rp(:) + r(:, ii)*sinj(ii)
5611 rpp(:) = rpp(:) + r(:, ii)*cosj(ii)
5614 nv = nv/sqrt(sum(nv**2))
5618 rpxpp = sqrt(sum(uv**2))
5623 dnvp(:, i) = uv - nv*sum(uv*nv)
5627 dnvpp(:, i) = uv - nv*sum(uv*nv)
5630 nforce(:, :, ii) = dnvp(:, :)*sinj(ii) + dnvpp(:, :)*cosj(ii)
5635 z(ii) = sum(r(:, ii)*nv(:))
5641 zforce(ii, jj, :) = nv
5643 zforce(ii, jj, :) = 0._dp
5647 zforce(ii, jj, i) = zforce(ii, jj, i) + r(j, ii)*nforce(j, i, jj)
5653 IF (colvar%ring_puckering_param%iq == 0)
THEN
5655 svar = sqrt(sum(z**2))
5659 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*z(jj)
5662 CALL put_derivative(colvar, ii, ftmp)
5665 m = abs(colvar%ring_puckering_param%iq)
5667 IF (mod(nring, 2) == 0 .AND. colvar%ring_puckering_param%iq == nring/2)
THEN
5671 IF (mod(ii, 2) == 0)
THEN
5677 svar = svar*sqrt(kr)
5681 IF (mod(jj, 2) == 0)
THEN
5682 ftmp(:) = ftmp(:) - zforce(jj, ii, :)*sqrt(kr)
5684 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*sqrt(kr)
5687 CALL put_derivative(colvar, ii, -ftmp)
5690 cpassert(m <= (nring - 1)/2)
5694 a = a + z(ii)*cos(
twopi*m*(ii - 1)*kr)
5695 b = b - z(ii)*sin(
twopi*m*(ii - 1)*kr)
5697 a = a*sqrt(2._dp*kr)
5698 b = b*sqrt(2._dp*kr)
5699 IF (colvar%ring_puckering_param%iq > 0)
THEN
5701 svar = sqrt(a*a + b*b)
5707 IF (at >
pi/2._dp)
THEN
5708 svar = 2.5_dp*
pi - at
5710 svar = 0.5_dp*
pi - at
5718 ds = da*cos(
twopi*m*(ii - 1)*kr)
5719 ds = ds - db*sin(
twopi*m*(ii - 1)*kr)
5720 ftmp(:) = ftmp(:) + ds*sqrt(2._dp*kr)*zforce(ii, jj, :)
5722 CALL put_derivative(colvar, jj, ftmp)
5729 DEALLOCATE (r, z, cosj, sinj, nforce, zforce)
5731 END SUBROUTINE ring_puckering_colvar
5753 RECURSIVE FUNCTION rec_eval_grid(iw1, ncol, f_vals, v_count, &
5754 gp, grid_sp, step_size, istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
RESULT(k)
5755 INTEGER :: iw1, ncol
5756 REAL(
dp),
DIMENSION(:, :),
POINTER :: f_vals
5758 REAL(
dp),
DIMENSION(:),
POINTER :: gp, grid_sp
5759 REAL(
dp) :: step_size
5760 INTEGER :: istart, iend
5761 REAL(
dp),
DIMENSION(:, :),
POINTER :: s1v
5762 REAL(
dp),
DIMENSION(:),
POINTER :: s1
5763 INTEGER,
DIMENSION(:, :),
POINTER :: p_bounds
5765 INTEGER :: ifunc, nconf, k
5767 INTEGER :: count1, i
5770 IF (v_count .LT. ncol)
THEN
5771 count1 = v_count + 1
5772 DO i = p_bounds(1, count1), p_bounds(2, count1)
5773 gp(count1) = real(i, kind=
dp)*grid_sp(count1)
5774 k = rec_eval_grid(iw1, ncol, f_vals, count1, gp, grid_sp, step_size, &
5775 istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
5777 ELSE IF (v_count == ncol .AND. ifunc == 1)
THEN
5779 s1v(1, i) = real(i, kind=
dp)*step_size*exp(-lambda*dot_product(gp(:) - f_vals(:, i), &
5780 gp(:) - f_vals(:, i)))
5781 s1v(2, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5784 s1(i) = accurate_sum(s1v(i, :))
5786 WRITE (iw1,
'(5F10.5)') gp(:), s1(1)/s1(2)/real(nconf - 1,
dp)
5787 ELSE IF (v_count == ncol .AND. ifunc == 2)
THEN
5789 s1v(1, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5791 s1(1) = accurate_sum(s1v(1, :))
5793 WRITE (iw1,
'(5F10.5)') gp(:), -lambda*log(s1(1))
5795 END FUNCTION rec_eval_grid
5808 SUBROUTINE read_frames(frame_section, para_env, nr_frames, r_ref, n_atoms)
5810 TYPE(section_vals_type),
POINTER :: frame_section
5811 TYPE(mp_para_env_type),
POINTER :: para_env
5812 INTEGER,
INTENT(IN) :: nr_frames
5813 REAL(
dp),
DIMENSION(:, :),
POINTER :: r_ref
5814 INTEGER,
INTENT(OUT) :: n_atoms
5816 CHARACTER(LEN=default_path_length) :: filename
5817 CHARACTER(LEN=default_string_length) :: dummy_char
5818 INTEGER :: i, j, natom
5819 LOGICAL :: explicit, my_end
5820 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rptr
5821 TYPE(section_vals_type),
POINTER :: coord_section
5833 ALLOCATE (r_ref(3*natom, nr_frames))
5836 cpassert(3*natom ==
SIZE(r_ref, 1))
5840 i_rep_val=j, r_vals=rptr)
5841 r_ref((j - 1)*3 + 1:(j - 1)*3 + 3, i) = rptr(1:3)
5845 TYPE(cp_parser_type) :: parser
5847 cpassert(trim(filename) /=
"")
5849 CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
5852 CALL parser_get_object(parser, natom)
5855 ALLOCATE (r_ref(3*natom, nr_frames))
5858 cpassert(3*natom ==
SIZE(r_ref, 1))
5864 CALL cp_abort(__location__, &
5865 "Number of lines in XYZ format not equal to the number of atoms."// &
5866 " Error in XYZ format for COORD_A (CV rmsd). Very probably the"// &
5867 " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
5868 READ (parser%input_line, *) dummy_char, rptr(1:3)
5879 END SUBROUTINE read_frames
5890 SUBROUTINE wc_colvar(colvar, cell, subsys, particles, qs_env)
5891 TYPE(colvar_type),
POINTER :: colvar
5892 TYPE(cell_type),
POINTER :: cell
5893 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5894 TYPE(particle_type),
DIMENSION(:), &
5895 OPTIONAL,
POINTER :: particles
5896 TYPE(qs_environment_type),
POINTER,
OPTIONAL :: qs_env
5898 INTEGER :: od, h, oa
5899 REAL(dp) :: rod(3), roa(3), rh(3), &
5900 x, y, s(3), xv(3), dmin, amin
5901 INTEGER :: idmin, iamin, i, j
5902 TYPE(particle_list_type),
POINTER :: particles_i
5903 TYPE(particle_type),
DIMENSION(:), &
5904 POINTER :: my_particles
5905 TYPE(wannier_centres_type),
DIMENSION(:),
POINTER :: wc
5906 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
5907 INTEGER :: nwca, nwcd
5910 NULLIFY (particles_i, wc)
5912 cpassert(colvar%type_id == wc_colvar_id)
5913 IF (
PRESENT(particles))
THEN
5914 my_particles => particles
5916 cpassert(
PRESENT(subsys))
5917 CALL cp_subsys_get(subsys, particles=particles_i)
5918 my_particles => particles_i%els
5920 CALL get_qs_env(qs_env, wanniercentres=wc)
5921 rcut = colvar%Wc%rcut
5922 od = colvar%Wc%ids(1)
5923 h = colvar%Wc%ids(2)
5924 oa = colvar%Wc%ids(3)
5925 CALL get_coordinates(colvar, od, rod, my_particles)
5926 CALL get_coordinates(colvar, h, rh, my_particles)
5927 CALL get_coordinates(colvar, oa, roa, my_particles)
5928 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
5929 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
5932 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
5933 x = distance(rod - wc(1)%centres(:, j))
5934 y = distance(roa - wc(1)%centres(:, j))
5946 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
5947 amin = distance(rh - wc(1)%centres(:, wcai(1)))
5952 x = distance(rh - wc(1)%centres(:, wcdi(i)))
5959 x = distance(rh - wc(1)%centres(:, wcai(i)))
5971 colvar%ss = wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
5981 REAL(dp) function distance(rij)
5982 REAL(dp),
INTENT(in) :: rij(3)
5984 s = matmul(cell%h_inv, rij)
5986 xv = matmul(cell%hmat, s)
5987 distance = sqrt(dot_product(xv, xv))
5988 END FUNCTION distance
5990 END SUBROUTINE wc_colvar
6001 SUBROUTINE hbp_colvar(colvar, cell, subsys, particles, qs_env)
6002 TYPE(colvar_type),
POINTER :: colvar
6003 TYPE(cell_type),
POINTER :: cell
6004 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
6005 TYPE(particle_type),
DIMENSION(:), &
6006 OPTIONAL,
POINTER :: particles
6007 TYPE(qs_environment_type),
POINTER,
OPTIONAL :: qs_env
6009 INTEGER :: od, h, oa
6010 REAL(dp) :: rod(3), roa(3), rh(3), &
6011 x, y, s(3), xv(3), dmin, amin
6012 INTEGER :: idmin, iamin, i, j, il, output_unit
6013 TYPE(particle_list_type),
POINTER :: particles_i
6014 TYPE(particle_type),
DIMENSION(:), &
6015 POINTER :: my_particles
6016 TYPE(wannier_centres_type), &
6017 DIMENSION(:),
POINTER :: wc
6018 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
6019 INTEGER :: nwca, nwcd
6022 NULLIFY (particles_i, wc)
6023 output_unit = cp_logger_get_default_io_unit()
6025 cpassert(colvar%type_id == hbp_colvar_id)
6026 IF (
PRESENT(particles))
THEN
6027 my_particles => particles
6029 cpassert(
PRESENT(subsys))
6030 CALL cp_subsys_get(subsys, particles=particles_i)
6031 my_particles => particles_i%els
6033 CALL get_qs_env(qs_env, wanniercentres=wc)
6034 rcut = colvar%HBP%rcut
6035 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
6036 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
6038 DO il = 1, colvar%HBP%nPoints
6039 od = colvar%HBP%ids(il, 1)
6040 h = colvar%HBP%ids(il, 2)
6041 oa = colvar%HBP%ids(il, 3)
6042 CALL get_coordinates(colvar, od, rod, my_particles)
6043 CALL get_coordinates(colvar, h, rh, my_particles)
6044 CALL get_coordinates(colvar, oa, roa, my_particles)
6047 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
6048 x = distance(rod - wc(1)%centres(:, j))
6049 y = distance(roa - wc(1)%centres(:, j))
6061 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
6062 amin = distance(rh - wc(1)%centres(:, wcai(1)))
6067 x = distance(rh - wc(1)%centres(:, wcdi(i)))
6074 x = distance(rh - wc(1)%centres(:, wcai(i)))
6080 colvar%HBP%ewc(il) = colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6081 colvar%ss = colvar%ss + colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6083 IF (output_unit > 0)
THEN
6084 DO il = 1, colvar%HBP%nPoints
6085 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP| = ", colvar%HBP%ewc(il)
6087 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP|\theta(x) = ", colvar%ss
6098 REAL(dp) function distance(rij)
6099 REAL(dp),
INTENT(in) :: rij(3)
6101 s = matmul(cell%h_inv, rij)
6103 xv = matmul(cell%hmat, s)
6104 distance = sqrt(dot_product(xv, xv))
6105 END FUNCTION distance
6107 END SUBROUTINE hbp_colvar
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Handles all functions related to the CELL.
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
subroutine, public colvar_eval_glob_f(icolvar, force_env)
evaluates the derivatives (dsdr) given and due to the given colvar
recursive subroutine, public colvar_read(colvar, icol, colvar_section, para_env)
reads a colvar from the input
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Initialize the collective variables types.
integer, parameter, public ring_puckering_colvar_id
integer, parameter, public population_colvar_id
integer, parameter, public do_clv_geo_center
integer, parameter, public distance_from_path_colvar_id
integer, parameter, public rmsd_colvar_id
integer, parameter, public mindist_colvar_id
integer, parameter, public wc_colvar_id
integer, parameter, public acid_hyd_dist_colvar_id
integer, parameter, public xyz_outerdiag_colvar_id
integer, parameter, public do_clv_xz
integer, parameter, public plane_plane_angle_colvar_id
subroutine, public colvar_create(colvar, colvar_id)
initializes a colvar_param type
integer, parameter, public plane_distance_colvar_id
integer, parameter, public combine_colvar_id
integer, parameter, public gyration_colvar_id
integer, parameter, public hbp_colvar_id
integer, parameter, public rotation_colvar_id
integer, parameter, public hydronium_dist_colvar_id
integer, parameter, public coord_colvar_id
integer, parameter, public do_clv_fix_point
integer, parameter, public do_clv_z
subroutine, public eval_point_pos(point, particles, r)
Evaluate the position of the geometrical point.
integer, parameter, public plane_def_atoms
integer, parameter, public do_clv_yz
integer, parameter, public dfunct_colvar_id
integer, parameter, public angle_colvar_id
integer, parameter, public qparm_colvar_id
subroutine, public eval_point_der(points, i, dsdr, f)
Evaluate the position of the geometrical point.
subroutine, public eval_point_mass(point, particles, m)
...
integer, parameter, public dist_colvar_id
subroutine, public colvar_setup(colvar)
Finalize the setup of the collective variable.
integer, parameter, public do_clv_xy
integer, parameter, public u_colvar_id
integer, parameter, public hydronium_shell_colvar_id
integer, parameter, public torsion_colvar_id
integer, parameter, public do_clv_y
integer, parameter, public plane_def_vec
integer, parameter, public xyz_diag_colvar_id
integer, parameter, public reaction_path_colvar_id
integer, parameter, public do_clv_x
integer, parameter, public acid_hyd_shell_colvar_id
subroutine, public check_fixed_atom_cns_colv(fixd_list, colvar)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Interface for the force calculations.
integer, parameter, public use_mixed_force
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
This public domain function parser module is intended for applications where a set of mathematical ex...
real(rn) function, public evalf(i, Val)
...
integer, public evalerrtype
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
subroutine, public finalizef()
...
subroutine, public initf(n)
...
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
integer, parameter, public maxfac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Collection of simple mathematical functions and subroutines.
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Utility routines for the memory handling.
Interface to the message passing library MPI.
subroutine, public get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, force_eval_embed)
performs mapping of the subsystems of different force_eval
Define the molecule kind structure types and the corresponding functionality.
represent a simple array based list of the given type
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Defines functions to perform rmsd in 3D.
subroutine, public rmsd3(particle_set, r, r0, output_unit, weights, my_val, rotate, transl, rot, drmsd3)
Computes the RMSD in 3D. Provides also derivatives.
Calculate spherical harmonics.
real(kind=dp) function, public legendre(x, l, m)
...
real(kind=dp) function, public dlegendre(x, l, m)
...
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
defines the type needed for computing wannier states expectations