99#include "./base/base_uses.f90"
104 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'colvar_methods'
105 REAL(KIND=
dp),
PRIVATE,
PARAMETER :: tolerance_acos = 1.0e-5_dp
124 RECURSIVE SUBROUTINE colvar_read(colvar, icol, colvar_section, para_env, cell)
126 INTEGER,
INTENT(IN) :: icol
129 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
131 CHARACTER(len=*),
PARAMETER :: routinen =
'colvar_read'
133 CHARACTER(LEN=3) :: fmid
134 CHARACTER(LEN=7) :: tag, tag_comp, tag_comp1, tag_comp2
135 CHARACTER(LEN=default_path_length) :: path_function
136 CHARACTER(LEN=default_string_length) :: tmpstr, tmpstr2
137 CHARACTER(LEN=default_string_length), &
138 DIMENSION(:),
POINTER :: c_kinds, my_par
139 INTEGER :: handle, i, iatm, icomponent, iend, &
140 ifunc, ii, isize, istart, iw, iw1, j, &
141 k, kk, n_var, n_var_k, ncol, ndim, &
143 INTEGER,
DIMENSION(:),
POINTER :: iatms
144 INTEGER,
DIMENSION(:, :),
POINTER :: p_bounds
145 LOGICAL :: check, use_mixed_energy
146 LOGICAL,
DIMENSION(26) :: my_subsection
147 REAL(
dp),
DIMENSION(:),
POINTER :: s1, wei, weights
148 REAL(
dp),
DIMENSION(:, :),
POINTER :: p_range, s1v
149 REAL(kind=
dp),
DIMENSION(1) :: my_val
150 REAL(kind=
dp),
DIMENSION(:),
POINTER :: g_range, grid_point, grid_sp, my_vals, &
156 TYPE(
section_vals_type),
POINTER :: acid_hyd_dist_section, acid_hyd_shell_section, &
157 angle_section, colvar_subsection, combine_section, coordination_section, dfunct_section, &
158 distance_from_path_section, distance_section, frame_section, gyration_section, &
159 hbp_section, hydronium_dist_section, hydronium_shell_section, mindist_section, &
160 path_section, plane_dist_section, plane_plane_angle_section, plane_sections, &
161 point_section, population_section, qparm_section, reaction_path_section, &
162 ring_puckering_section, rmsd_section, rotation_section, torsion_section, u_section, &
163 wc_section, wrk_section
166 CALL timeset(routinen, handle)
167 NULLIFY (logger, c_kinds, iatms)
169 my_subsection = .false.
177 plane_plane_angle_section &
184 acid_hyd_shell_section &
187 can_return_null=.true.)
188 distance_from_path_section &
190 i_rep_section=icol, can_return_null=.true.)
192 can_return_null=.true.)
201 ring_puckering_section &
215 IF (
ASSOCIATED(reaction_path_section))
THEN
217 explicit=my_subsection(10))
219 IF (
ASSOCIATED(distance_from_path_section))
THEN
221 explicit=my_subsection(16))
223 IF (
ASSOCIATED(combine_section))
THEN
228 explicit=my_subsection(13))
237 explicit=my_subsection(22))
244 cpassert(count(my_subsection) == 1)
245 cpassert(.NOT.
ASSOCIATED(colvar))
247 IF (my_subsection(1))
THEN
249 wrk_section => distance_section
251 CALL colvar_check_points(colvar, distance_section, cell)
253 colvar%dist_param%i_at = iatms(1)
254 colvar%dist_param%j_at = iatms(2)
257 ELSE IF (my_subsection(2))
THEN
259 wrk_section => angle_section
261 CALL colvar_check_points(colvar, angle_section, cell)
263 colvar%angle_param%i_at_angle = iatms
264 ELSE IF (my_subsection(3))
THEN
266 wrk_section => torsion_section
268 CALL colvar_check_points(colvar, torsion_section, cell)
270 colvar%torsion_param%i_at_tors = iatms
271 colvar%torsion_param%o0 = 0.0_dp
272 ELSE IF (my_subsection(4))
THEN
274 wrk_section => coordination_section
276 CALL colvar_check_points(colvar, coordination_section, cell)
277 NULLIFY (colvar%coord_param%i_at_from, colvar%coord_param%c_kinds_from)
278 NULLIFY (colvar%coord_param%i_at_to, colvar%coord_param%c_kinds_to)
279 NULLIFY (colvar%coord_param%i_at_to_b, colvar%coord_param%c_kinds_to_b)
287 CALL reallocate(colvar%coord_param%i_at_from, 1, ndim +
SIZE(iatms))
288 colvar%coord_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
289 ndim = ndim +
SIZE(iatms)
291 colvar%coord_param%n_atoms_from = ndim
292 colvar%coord_param%use_kinds_from = .false.
299 CALL reallocate(colvar%coord_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
300 colvar%coord_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
301 ndim = ndim +
SIZE(c_kinds)
303 colvar%coord_param%n_atoms_from = 0
304 colvar%coord_param%use_kinds_from = .true.
307 CALL uppercase(colvar%coord_param%c_kinds_from(k))
317 CALL reallocate(colvar%coord_param%i_at_to, 1, ndim +
SIZE(iatms))
318 colvar%coord_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
319 ndim = ndim +
SIZE(iatms)
321 colvar%coord_param%n_atoms_to = ndim
322 colvar%coord_param%use_kinds_to = .false.
329 CALL reallocate(colvar%coord_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
330 colvar%coord_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
331 ndim = ndim +
SIZE(c_kinds)
333 colvar%coord_param%n_atoms_to = 0
334 colvar%coord_param%use_kinds_to = .true.
337 CALL uppercase(colvar%coord_param%c_kinds_to(k))
348 IF (n_var /= 0 .OR. n_var_k /= 0)
THEN
349 colvar%coord_param%do_chain = .true.
354 CALL reallocate(colvar%coord_param%i_at_to_b, 1, ndim +
SIZE(iatms))
355 colvar%coord_param%i_at_to_b(ndim + 1:ndim +
SIZE(iatms)) = iatms
356 ndim = ndim +
SIZE(iatms)
358 colvar%coord_param%n_atoms_to_b = ndim
359 colvar%coord_param%use_kinds_to_b = .false.
363 cpassert(n_var_k > 0)
366 CALL reallocate(colvar%coord_param%c_kinds_to_b, 1, ndim +
SIZE(c_kinds))
367 colvar%coord_param%c_kinds_to_b(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
368 ndim = ndim +
SIZE(c_kinds)
370 colvar%coord_param%n_atoms_to_b = 0
371 colvar%coord_param%use_kinds_to_b = .true.
374 CALL uppercase(colvar%coord_param%c_kinds_to_b(k))
382 colvar%coord_param%do_chain = .false.
383 colvar%coord_param%n_atoms_to_b = 0
384 colvar%coord_param%use_kinds_to_b = .false.
385 NULLIFY (colvar%coord_param%i_at_to_b)
386 NULLIFY (colvar%coord_param%c_kinds_to_b)
387 colvar%coord_param%nncrd_b = 0
388 colvar%coord_param%ndcrd_b = 0
389 colvar%coord_param%r_0_b = 0._dp
392 ELSE IF (my_subsection(5))
THEN
394 wrk_section => plane_dist_section
396 CALL colvar_check_points(colvar, plane_dist_section, cell)
398 cpassert(
SIZE(iatms) == 3)
399 colvar%plane_distance_param%plane = iatms
401 colvar%plane_distance_param%point = iatm
403 ELSE IF (my_subsection(6))
THEN
405 wrk_section => rotation_section
407 CALL colvar_check_points(colvar, rotation_section, cell)
408 CALL section_vals_val_get(rotation_section,
"P1_BOND1", i_val=colvar%rotation_param%i_at1_bond1)
409 CALL section_vals_val_get(rotation_section,
"P2_BOND1", i_val=colvar%rotation_param%i_at2_bond1)
410 CALL section_vals_val_get(rotation_section,
"P1_BOND2", i_val=colvar%rotation_param%i_at1_bond2)
411 CALL section_vals_val_get(rotation_section,
"P2_BOND2", i_val=colvar%rotation_param%i_at2_bond2)
412 ELSE IF (my_subsection(7))
THEN
414 wrk_section => dfunct_section
416 CALL colvar_check_points(colvar, dfunct_section, cell)
418 colvar%dfunct_param%i_at_dfunct = iatms
421 ELSE IF (my_subsection(8))
THEN
423 wrk_section => qparm_section
425 CALL colvar_check_points(colvar, qparm_section, cell)
428 CALL section_vals_val_get(qparm_section,
"INCLUDE_IMAGES", l_val=colvar%qparm_param%include_images)
431 NULLIFY (colvar%qparm_param%i_at_from)
432 NULLIFY (colvar%qparm_param%i_at_to)
437 CALL reallocate(colvar%qparm_param%i_at_from, 1, ndim +
SIZE(iatms))
438 colvar%qparm_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
439 ndim = ndim +
SIZE(iatms)
441 colvar%qparm_param%n_atoms_from = ndim
447 CALL reallocate(colvar%qparm_param%i_at_to, 1, ndim +
SIZE(iatms))
448 colvar%qparm_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
449 ndim = ndim +
SIZE(iatms)
451 colvar%qparm_param%n_atoms_to = ndim
452 ELSE IF (my_subsection(9))
THEN
455 NULLIFY (colvar%hydronium_shell_param%i_oxygens)
456 NULLIFY (colvar%hydronium_shell_param%i_hydrogens)
458 colvar%hydronium_shell_param%n_oxygens, &
459 colvar%hydronium_shell_param%n_hydrogens, &
460 colvar%hydronium_shell_param%i_oxygens, &
461 colvar%hydronium_shell_param%i_hydrogens)
462 ELSE IF (my_subsection(10) .OR. my_subsection(16))
THEN
464 IF (my_subsection(10))
THEN
465 path_section => reaction_path_section
469 ELSE IF (my_subsection(16))
THEN
470 path_section => distance_from_path_section
475 colvar%use_points = .false.
477 CALL section_vals_val_get(path_section,
"DISTANCES_RMSD", l_val=colvar%reaction_path_param%dist_rmsd)
479 IF (colvar%reaction_path_param%dist_rmsd .AND. colvar%reaction_path_param%rmsd)
THEN
480 cpabort(
"CV REACTION PATH: only one between DISTANCES_RMSD and RMSD can be used ")
482 IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd)
THEN
483 NULLIFY (colvar%reaction_path_param%i_rmsd, colvar%reaction_path_param%r_ref)
487 colvar%reaction_path_param%nr_frames = nr_frame
488 CALL read_frames(frame_section, para_env, nr_frame, colvar%reaction_path_param%r_ref, &
489 colvar%reaction_path_param%n_components)
491 IF (colvar%reaction_path_param%subset ==
rmsd_all)
THEN
492 ALLOCATE (colvar%reaction_path_param%i_rmsd(colvar%reaction_path_param%n_components))
493 DO i = 1, colvar%reaction_path_param%n_components
494 colvar%reaction_path_param%i_rmsd(i) = i
496 ELSE IF (colvar%reaction_path_param%subset ==
rmsd_list)
THEN
504 CALL reallocate(colvar%reaction_path_param%i_rmsd, 1, ndim +
SIZE(iatms))
505 colvar%reaction_path_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
506 ndim = ndim +
SIZE(iatms)
508 colvar%reaction_path_param%n_components = ndim
510 cpabort(
"CV REACTION PATH: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
514 CALL section_vals_val_get(path_section,
"ALIGN_FRAMES", l_val=colvar%reaction_path_param%align_frames)
518 ALLOCATE (colvar%reaction_path_param%colvar_p(ncol))
521 NULLIFY (colvar%reaction_path_param%colvar_p(i)%colvar)
522 CALL colvar_read(colvar%reaction_path_param%colvar_p(i)%colvar, i, colvar_subsection, para_env, cell)
525 cpabort(
"CV REACTION PATH: the number of CV to define the path must be >0 ")
527 colvar%reaction_path_param%n_components = ncol
530 CALL section_vals_val_get(path_section,
"STEP_SIZE", r_val=colvar%reaction_path_param%step_size)
531 iend = ceiling(max(range(1), range(2))/colvar%reaction_path_param%step_size)
532 istart = floor(min(range(1), range(2))/colvar%reaction_path_param%step_size)
533 colvar%reaction_path_param%function_bounds(1) = istart
534 colvar%reaction_path_param%function_bounds(2) = iend
535 colvar%reaction_path_param%nr_frames = 2
536 ALLOCATE (colvar%reaction_path_param%f_vals(ncol, istart:iend))
539 check = (ncol ==
SIZE(colvar%reaction_path_param%colvar_p))
544 CALL compress(path_function, full=.true.)
545 CALL parsef(i, trim(path_function), my_par)
547 my_val = real(j, kind=
dp)*colvar%reaction_path_param%step_size
548 colvar%reaction_path_param%f_vals(i, j) =
evalf(i, my_val)
554 "MAP", middle_name=fmid, extension=
".dat", file_status=
"REPLACE")
557 ALLOCATE (grid_sp(ncol))
562 cpassert(ncol ==
SIZE(grid_sp))
563 ALLOCATE (p_range(2, ncol))
564 ALLOCATE (p_bounds(2, ncol))
567 p_range(:, i) = g_range(:)
568 p_bounds(2, i) = ceiling(max(p_range(1, i), p_range(2, i))/grid_sp(i))
569 p_bounds(1, i) = floor(min(p_range(1, i), p_range(2, i))/grid_sp(i))
571 ALLOCATE (s1v(2, istart:iend))
573 ALLOCATE (grid_point(ncol))
575 kk = rec_eval_grid(iw1, ncol, colvar%reaction_path_param%f_vals, v_count, &
576 grid_point, grid_sp, colvar%reaction_path_param%step_size, istart, &
577 iend, s1v, s1, p_bounds, colvar%reaction_path_param%lambda, ifunc=ifunc, &
578 nconf=colvar%reaction_path_param%nr_frames)
581 DEALLOCATE (p_bounds)
584 DEALLOCATE (grid_point)
590 ELSE IF (my_subsection(11))
THEN
593 colvar%use_points = .false.
596 ALLOCATE (colvar%combine_cvs_param%colvar_p(ncol))
599 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
601 WRITE (iw,
'( A )')
' '// &
602 '**********************************************************************'
603 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
604 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| COMBINATION OF THE FOLLOWING COLVARS:'
607 "PRINT%PROGRAM_RUN_INFO")
610 NULLIFY (colvar%combine_cvs_param%colvar_p(i)%colvar)
611 CALL colvar_read(colvar%combine_cvs_param%colvar_p(i)%colvar, i, colvar_subsection, para_env, cell)
615 CALL compress(colvar%combine_cvs_param%function, full=.true.)
618 ALLOCATE (colvar%combine_cvs_param%variables(
SIZE(my_par)))
619 colvar%combine_cvs_param%variables = my_par
621 IF (
SIZE(my_par) /= ncol)
THEN
622 CALL cp_abort(__location__, &
623 "Number of defined COLVAR for COMBINE_COLVAR is different from the "// &
624 "number of variables! It is not possible to define COLVARs in a COMBINE_COLVAR "// &
625 "and avoid their usage in the combininig function!")
628 ALLOCATE (colvar%combine_cvs_param%c_parameters(0))
631 isize =
SIZE(colvar%combine_cvs_param%c_parameters)
633 CALL reallocate(colvar%combine_cvs_param%c_parameters, 1, isize +
SIZE(my_par))
634 colvar%combine_cvs_param%c_parameters(isize + 1:isize +
SIZE(my_par)) = my_par
636 ALLOCATE (colvar%combine_cvs_param%v_parameters(0))
639 isize =
SIZE(colvar%combine_cvs_param%v_parameters)
641 CALL reallocate(colvar%combine_cvs_param%v_parameters, 1, isize +
SIZE(my_vals))
642 colvar%combine_cvs_param%v_parameters(isize + 1:isize +
SIZE(my_vals)) = my_vals
647 ELSE IF (my_subsection(12))
THEN
649 wrk_section => population_section
651 CALL colvar_check_points(colvar, population_section, cell)
653 NULLIFY (colvar%population_param%i_at_from, colvar%population_param%c_kinds_from)
654 NULLIFY (colvar%population_param%i_at_to, colvar%population_param%c_kinds_to)
663 CALL reallocate(colvar%population_param%i_at_from, 1, ndim +
SIZE(iatms))
664 colvar%population_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
665 ndim = ndim +
SIZE(iatms)
667 colvar%population_param%n_atoms_from = ndim
668 colvar%population_param%use_kinds_from = .false.
675 CALL reallocate(colvar%population_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
676 colvar%population_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
677 ndim = ndim +
SIZE(c_kinds)
679 colvar%population_param%n_atoms_from = 0
680 colvar%population_param%use_kinds_from = .true.
683 CALL uppercase(colvar%population_param%c_kinds_from(k))
693 CALL reallocate(colvar%population_param%i_at_to, 1, ndim +
SIZE(iatms))
694 colvar%population_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
695 ndim = ndim +
SIZE(iatms)
697 colvar%population_param%n_atoms_to = ndim
698 colvar%population_param%use_kinds_to = .false.
705 CALL reallocate(colvar%population_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
706 colvar%population_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
707 ndim = ndim +
SIZE(c_kinds)
709 colvar%population_param%n_atoms_to = 0
710 colvar%population_param%use_kinds_to = .true.
713 CALL uppercase(colvar%population_param%c_kinds_to(k))
722 ELSE IF (my_subsection(13))
THEN
724 wrk_section => plane_plane_angle_section
726 CALL colvar_check_points(colvar, plane_plane_angle_section, cell)
731 cpabort(
"PLANE_PLANE_ANGLE Colvar section: Two PLANE sections must be provided!")
735 i_val=colvar%plane_plane_angle_param%plane1%type_of_def)
736 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_vec)
THEN
739 colvar%plane_plane_angle_param%plane1%normal_vec = s1
740 IF (
PRESENT(cell))
THEN
741 IF (
ASSOCIATED(cell))
THEN
748 colvar%plane_plane_angle_param%plane1%points = iatms
753 i_val=colvar%plane_plane_angle_param%plane2%type_of_def)
754 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_vec)
THEN
757 colvar%plane_plane_angle_param%plane2%normal_vec = s1
758 IF (
PRESENT(cell))
THEN
759 IF (
ASSOCIATED(cell))
THEN
766 colvar%plane_plane_angle_param%plane2%points = iatms
768 ELSE IF (my_subsection(14))
THEN
770 wrk_section => gyration_section
772 CALL colvar_check_points(colvar, gyration_section, cell)
774 NULLIFY (colvar%gyration_param%i_at, colvar%gyration_param%c_kinds)
783 CALL reallocate(colvar%gyration_param%i_at, 1, ndim +
SIZE(iatms))
784 colvar%gyration_param%i_at(ndim + 1:ndim +
SIZE(iatms)) = iatms
785 ndim = ndim +
SIZE(iatms)
787 colvar%gyration_param%n_atoms = ndim
788 colvar%gyration_param%use_kinds = .false.
795 CALL reallocate(colvar%gyration_param%c_kinds, 1, ndim +
SIZE(c_kinds))
796 colvar%gyration_param%c_kinds(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
797 ndim = ndim +
SIZE(c_kinds)
799 colvar%gyration_param%n_atoms = 0
800 colvar%gyration_param%use_kinds = .true.
803 CALL uppercase(colvar%gyration_param%c_kinds(k))
806 ELSE IF (my_subsection(15))
THEN
808 wrk_section => rmsd_section
811 NULLIFY (colvar%rmsd_param%i_rmsd, colvar%rmsd_param%r_ref, colvar%rmsd_param%weights)
816 colvar%rmsd_param%nr_frames = nr_frame
818 cpassert(nr_frame >= 1 .AND. nr_frame <= 2)
819 CALL read_frames(frame_section, para_env, nr_frame, colvar%rmsd_param%r_ref, &
820 colvar%rmsd_param%n_atoms)
821 ALLOCATE (colvar%rmsd_param%weights(colvar%rmsd_param%n_atoms))
822 colvar%rmsd_param%weights = 0.0_dp
824 IF (colvar%rmsd_param%subset ==
rmsd_all)
THEN
825 ALLOCATE (colvar%rmsd_param%i_rmsd(colvar%rmsd_param%n_atoms))
826 DO i = 1, colvar%rmsd_param%n_atoms
827 colvar%rmsd_param%i_rmsd(i) = i
829 ELSE IF (colvar%rmsd_param%subset ==
rmsd_list)
THEN
837 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
838 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
839 ndim = ndim +
SIZE(iatms)
841 colvar%rmsd_param%n_atoms = ndim
843 cpabort(
"CV RMSD: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
852 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
853 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
854 ndim = ndim +
SIZE(iatms)
856 colvar%rmsd_param%n_atoms = ndim
858 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of atoms needs to be provided ")
867 weights(ndim + 1:ndim +
SIZE(wei)) = wei
868 ndim = ndim +
SIZE(wei)
870 IF (ndim /= colvar%rmsd_param%n_atoms)
THEN
871 CALL cp_abort(__location__,
"CV RMSD: list of atoms and list of "// &
872 "weights need to contain same number of entries. ")
875 ii = colvar%rmsd_param%i_rmsd(i)
876 colvar%rmsd_param%weights(ii) = weights(i)
880 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of weights need to be provided. ")
884 cpabort(
"CV RMSD: unknown SUBSET_TYPE.")
889 ELSE IF (my_subsection(17))
THEN
891 wrk_section => xyz_diag_section
893 CALL colvar_check_points(colvar, wrk_section, cell)
897 CALL section_vals_val_get(wrk_section,
"ABSOLUTE_POSITION", l_val=colvar%xyz_diag_param%use_absolute_position)
898 colvar%xyz_diag_param%i_atom = iatm
899 colvar%xyz_diag_param%component = icomponent
900 ELSE IF (my_subsection(18))
THEN
902 wrk_section => xyz_outerdiag_section
904 CALL colvar_check_points(colvar, wrk_section, cell)
906 colvar%xyz_outerdiag_param%i_atoms = iatms
908 colvar%xyz_outerdiag_param%components(1) = icomponent
910 colvar%xyz_outerdiag_param%components(2) = icomponent
912 ELSE IF (my_subsection(19))
THEN
914 wrk_section => u_section
917 CALL section_vals_get(colvar%u_param%mixed_energy_section, explicit=use_mixed_energy)
918 IF (.NOT. use_mixed_energy)
NULLIFY (colvar%u_param%mixed_energy_section)
919 ELSE IF (my_subsection(20))
THEN
921 wrk_section => wc_section
923 CALL colvar_check_points(colvar, wc_section, cell)
927 colvar%Wc%ids = iatms
928 ELSE IF (my_subsection(21))
THEN
930 wrk_section => hbp_section
932 CALL colvar_check_points(colvar, hbp_section, cell)
938 ALLOCATE (colvar%HBP%ids(colvar%HBP%nPoints, 3))
939 ALLOCATE (colvar%HBP%ewc(colvar%HBP%nPoints))
940 DO i = 1, colvar%HBP%nPoints
942 colvar%HBP%ids(i, :) = iatms
944 ELSE IF (my_subsection(22))
THEN
948 colvar%ring_puckering_param%nring =
SIZE(iatms)
949 ALLOCATE (colvar%ring_puckering_param%atoms(
SIZE(iatms)))
950 colvar%ring_puckering_param%atoms = iatms
952 i_val=colvar%ring_puckering_param%iq)
954 ndim = colvar%ring_puckering_param%nring
956 cpabort(
"CV Ring Puckering: Ring size has to be 4 or larger. ")
958 ii = colvar%ring_puckering_param%iq
959 IF (abs(ii) == 1 .OR. ii < -(ndim - 1)/2 .OR. ii > ndim/2)
THEN
960 cpabort(
"CV Ring Puckering: Invalid coordinate number.")
962 ELSE IF (my_subsection(23))
THEN
964 wrk_section => mindist_section
966 CALL colvar_check_points(colvar, mindist_section, cell)
967 NULLIFY (colvar%mindist_param%i_dist_from, colvar%mindist_param%i_coord_from, &
968 colvar%mindist_param%k_coord_from, colvar%mindist_param%i_coord_to, &
969 colvar%mindist_param%k_coord_to)
971 colvar%mindist_param%n_dist_from =
SIZE(iatms)
972 ALLOCATE (colvar%mindist_param%i_dist_from(
SIZE(iatms)))
973 colvar%mindist_param%i_dist_from = iatms
980 CALL reallocate(colvar%mindist_param%i_coord_from, 1, ndim +
SIZE(iatms))
981 colvar%mindist_param%i_coord_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
982 ndim = ndim +
SIZE(iatms)
984 colvar%mindist_param%n_coord_from = ndim
985 colvar%mindist_param%use_kinds_from = .false.
992 CALL reallocate(colvar%mindist_param%k_coord_from, 1, ndim +
SIZE(c_kinds))
993 colvar%mindist_param%k_coord_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
994 ndim = ndim +
SIZE(c_kinds)
996 colvar%mindist_param%n_coord_from = 0
997 colvar%mindist_param%use_kinds_from = .true.
1000 CALL uppercase(colvar%mindist_param%k_coord_from(k))
1006 IF (n_var /= 0)
THEN
1010 CALL reallocate(colvar%mindist_param%i_coord_to, 1, ndim +
SIZE(iatms))
1011 colvar%mindist_param%i_coord_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
1012 ndim = ndim +
SIZE(iatms)
1014 colvar%mindist_param%n_coord_to = ndim
1015 colvar%mindist_param%use_kinds_to = .false.
1022 CALL reallocate(colvar%mindist_param%k_coord_to, 1, ndim +
SIZE(c_kinds))
1023 colvar%mindist_param%k_coord_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
1024 ndim = ndim +
SIZE(c_kinds)
1026 colvar%mindist_param%n_coord_to = 0
1027 colvar%mindist_param%use_kinds_to = .true.
1030 CALL uppercase(colvar%mindist_param%k_coord_to(k))
1039 ELSE IF (my_subsection(24))
THEN
1042 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_water)
1043 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_acid)
1044 NULLIFY (colvar%acid_hyd_dist_param%i_hydrogens)
1046 colvar%acid_hyd_dist_param%n_oxygens_water, &
1047 colvar%acid_hyd_dist_param%n_oxygens_acid, &
1048 colvar%acid_hyd_dist_param%n_hydrogens, &
1049 colvar%acid_hyd_dist_param%i_oxygens_water, &
1050 colvar%acid_hyd_dist_param%i_oxygens_acid, &
1051 colvar%acid_hyd_dist_param%i_hydrogens)
1052 ELSE IF (my_subsection(25))
THEN
1055 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_water)
1056 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_acid)
1057 NULLIFY (colvar%acid_hyd_shell_param%i_hydrogens)
1059 colvar%acid_hyd_shell_param%n_oxygens_water, &
1060 colvar%acid_hyd_shell_param%n_oxygens_acid, &
1061 colvar%acid_hyd_shell_param%n_hydrogens, &
1062 colvar%acid_hyd_shell_param%i_oxygens_water, &
1063 colvar%acid_hyd_shell_param%i_oxygens_acid, &
1064 colvar%acid_hyd_shell_param%i_hydrogens)
1065 ELSE IF (my_subsection(26))
THEN
1068 NULLIFY (colvar%hydronium_dist_param%i_oxygens)
1069 NULLIFY (colvar%hydronium_dist_param%i_hydrogens)
1071 colvar%hydronium_dist_param%n_oxygens, &
1072 colvar%hydronium_dist_param%n_hydrogens, &
1073 colvar%hydronium_dist_param%i_oxygens, &
1074 colvar%hydronium_dist_param%i_hydrogens)
1079 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
1082 IF (colvar%use_points) tag =
"POINTS:"
1085 WRITE (iw,
'( A )')
' '// &
1086 '----------------------------------------------------------------------'
1087 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
1090 SELECT CASE (colvar%type_id)
1092 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| ANGLE >>> '//tag, &
1093 colvar%angle_param%i_at_angle
1095 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| DISTANCE DIFFERENCE >>> '//tag, &
1096 colvar%dfunct_param%i_at_dfunct
1098 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE DISTANCE - PLANE >>> '//tag, &
1099 colvar%plane_distance_param%plane
1100 WRITE (iw,
'( A,T73,1I8)')
' COLVARS| PLANE DISTANCE - POINT >>> '//tag, &
1101 colvar%plane_distance_param%point
1103 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
1104 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1105 colvar%plane_plane_angle_param%plane1%points
1107 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1108 colvar%plane_plane_angle_param%plane1%normal_vec
1111 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
1112 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1113 colvar%plane_plane_angle_param%plane2%points
1115 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1116 colvar%plane_plane_angle_param%plane2%normal_vec
1119 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| TORSION >>> '//tag, &
1120 colvar%torsion_param%i_at_tors
1122 WRITE (iw,
'( A,T65,2I8)')
' COLVARS| BOND >>> '//tag, &
1123 colvar%dist_param%i_at, colvar%dist_param%j_at
1125 IF (colvar%coord_param%do_chain)
THEN
1126 WRITE (iw,
'( A)')
' COLVARS| COORDINATION CHAIN FC(from->to)*FC(to->to_B)>> '
1128 IF (colvar%coord_param%use_kinds_from)
THEN
1129 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> FROM KINDS', &
1130 adjustr(colvar%coord_param%c_kinds_from(kk) (1:10)), &
1131 kk=1,
SIZE(colvar%coord_param%c_kinds_from))
1133 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> FROM '//tag, &
1134 colvar%coord_param%i_at_from(kk), &
1135 kk=1,
SIZE(colvar%coord_param%i_at_from))
1137 IF (colvar%coord_param%use_kinds_to)
THEN
1138 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS', &
1139 adjustr(colvar%coord_param%c_kinds_to(kk) (1:10)), &
1140 kk=1,
SIZE(colvar%coord_param%c_kinds_to))
1142 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag, &
1143 colvar%coord_param%i_at_to(kk), &
1144 kk=1,
SIZE(colvar%coord_param%i_at_to))
1146 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%coord_param%r_0
1147 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%coord_param%nncrd
1148 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%coord_param%ndcrd
1149 IF (colvar%coord_param%do_chain)
THEN
1150 IF (colvar%coord_param%use_kinds_to_b)
THEN
1151 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS B', &
1152 adjustr(colvar%coord_param%c_kinds_to_b(kk) (1:10)), &
1153 kk=1,
SIZE(colvar%coord_param%c_kinds_to_b))
1155 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag//
' B', &
1156 colvar%coord_param%i_at_to_b(kk), &
1157 kk=1,
SIZE(colvar%coord_param%i_at_to_b))
1159 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0 B', colvar%coord_param%r_0_b
1160 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN B', colvar%coord_param%nncrd_b
1161 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND B', colvar%coord_param%ndcrd_b
1164 IF (colvar%population_param%use_kinds_from)
THEN
1165 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> FROM KINDS', &
1166 adjustr(colvar%population_param%c_kinds_from(kk) (1:10)), &
1167 kk=1,
SIZE(colvar%population_param%c_kinds_from))
1169 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> FROM '//tag, &
1170 colvar%population_param%i_at_from(kk), &
1171 kk=1,
SIZE(colvar%population_param%i_at_from))
1173 IF (colvar%population_param%use_kinds_to)
THEN
1174 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> TO KINDS', &
1175 adjustr(colvar%population_param%c_kinds_to(kk) (1:10)), &
1176 kk=1,
SIZE(colvar%population_param%c_kinds_to))
1178 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> TO '//tag, &
1179 colvar%population_param%i_at_to(kk), &
1180 kk=1,
SIZE(colvar%population_param%i_at_to))
1182 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%population_param%r_0
1183 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%population_param%nncrd
1184 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%population_param%ndcrd
1185 WRITE (iw,
'( A,T71,I10)')
' COLVARS| N0', colvar%population_param%n0
1186 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| SIGMA', colvar%population_param%sigma
1188 IF (colvar%gyration_param%use_kinds)
THEN
1189 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| Gyration Radius >>> KINDS', &
1190 adjustr(colvar%gyration_param%c_kinds(kk) (1:10)), &
1191 kk=1,
SIZE(colvar%gyration_param%c_kinds))
1193 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Gyration Radius >>> ATOMS '//tag, &
1194 colvar%gyration_param%i_at(kk), &
1195 kk=1,
SIZE(colvar%gyration_param%i_at))
1198 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 1 >>> '//tag, &
1199 colvar%rotation_param%i_at1_bond1
1200 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 1 >>> '//tag, &
1201 colvar%rotation_param%i_at2_bond1
1202 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 2 >>> '//tag, &
1203 colvar%rotation_param%i_at1_bond2
1204 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 2 >>> '//tag, &
1205 colvar%rotation_param%i_at2_bond2
1207 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> FROM '//tag, &
1208 colvar%qparm_param%i_at_from(kk), &
1209 kk=1,
SIZE(colvar%qparm_param%i_at_from))
1210 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> TO '//tag, &
1211 colvar%qparm_param%i_at_to(kk), &
1212 kk=1,
SIZE(colvar%qparm_param%i_at_to))
1213 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RCUT', colvar%qparm_param%rcut
1214 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RSTART', colvar%qparm_param%rstart
1215 WRITE (iw,
'( A,T71,L10)')
' COLVARS| INCLUDE IMAGES', colvar%qparm_param%include_images
1217 WRITE (iw,
'( A,T71,I10)')
' COLVARS| L', colvar%qparm_param%l
1219 WRITE (iw,
'( A)')
' COLVARS| COMBINING FUNCTION : '// &
1220 trim(colvar%combine_cvs_param%function)
1221 WRITE (iw,
'( A)', advance=
"NO")
' COLVARS| VARIABLES : '
1222 DO i = 1,
SIZE(colvar%combine_cvs_param%variables)
1223 WRITE (iw,
'( A)', advance=
"NO") &
1224 trim(colvar%combine_cvs_param%variables(i))//
" "
1227 WRITE (iw,
'( A)')
' COLVARS| DEFINED PARAMETERS [label] [value]:'
1228 DO i = 1,
SIZE(colvar%combine_cvs_param%c_parameters)
1229 WRITE (iw,
'( A,A7,F9.3)')
' ', &
1230 trim(colvar%combine_cvs_param%c_parameters(i)), colvar%combine_cvs_param%v_parameters(i)
1232 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| ERROR ON DERIVATIVE EVALUATION', &
1233 colvar%combine_cvs_param%lerr
1234 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| DX', &
1235 colvar%combine_cvs_param%dx
1237 cpwarn(
"Description header for REACTION_PATH COLVAR missing!")
1239 cpwarn(
"Description header for REACTION_PATH COLVAR missing!")
1241 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_shell_param%poh
1242 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_shell_param%qoh
1243 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%hydronium_shell_param%poo
1244 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%hydronium_shell_param%qoo
1245 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%hydronium_shell_param%roo
1246 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_shell_param%roh
1247 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_shell_param%nh
1248 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%hydronium_shell_param%lambda
1250 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_dist_param%poh
1251 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_dist_param%qoh
1252 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_dist_param%roh
1253 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%hydronium_dist_param%pm
1254 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%hydronium_dist_param%qm
1255 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_dist_param%nh
1256 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PF', colvar%hydronium_dist_param%pf
1257 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QF', colvar%hydronium_dist_param%qf
1258 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NN', colvar%hydronium_dist_param%nn
1260 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_dist_param%paoh
1261 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_dist_param%qaoh
1262 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_dist_param%pwoh
1263 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_dist_param%qwoh
1264 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_dist_param%pcut
1265 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_dist_param%qcut
1266 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_dist_param%raoh
1267 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_dist_param%rwoh
1268 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_dist_param%nc
1269 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_dist_param%lambda
1271 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_shell_param%paoh
1272 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_shell_param%qaoh
1273 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_shell_param%pwoh
1274 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_shell_param%qwoh
1275 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%acid_hyd_shell_param%poo
1276 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%acid_hyd_shell_param%qoo
1277 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%acid_hyd_shell_param%pm
1278 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%acid_hyd_shell_param%qm
1279 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_shell_param%pcut
1280 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_shell_param%qcut
1281 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_shell_param%raoh
1282 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_shell_param%rwoh
1283 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%acid_hyd_shell_param%roo
1284 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%acid_hyd_shell_param%nh
1285 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_shell_param%nc
1286 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_shell_param%lambda
1288 cpwarn(
"Description header for RMSD COLVAR missing!")
1290 NULLIFY (section, keyword, enum)
1294 tag_comp = trim(
enum_i2c(enum, colvar%xyz_diag_param%component))
1297 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| POSITION ('//trim(tag_comp) &
1298 //
') >>> '//tag, colvar%xyz_diag_param%i_atom
1300 NULLIFY (section, keyword, enum)
1304 tag_comp1 = trim(
enum_i2c(enum, colvar%xyz_outerdiag_param%components(1)))
1307 tag_comp2 = trim(
enum_i2c(enum, colvar%xyz_outerdiag_param%components(2)))
1310 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| CROSS TERM POSITION ('//trim(tag_comp1) &
1311 //
" * "//trim(tag_comp2)//
') >>> '//tag, colvar%xyz_outerdiag_param%i_atoms
1313 WRITE (iw,
'( A,T77,A4)')
' COLVARS| ENERGY >>> '//tag,
'all!'
1315 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| Wc >>> RCUT: ', &
1317 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| Wc >>> '//tag, &
1320 WRITE (iw,
'( A,T57,I8)')
' COLVARS| HBP >>> NPOINTS', &
1322 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1324 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1326 DO i = 1, colvar%HBP%nPoints
1327 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| HBP >>> '//tag, &
1328 colvar%HBP%ids(i, :)
1331 WRITE (iw,
'( A,T57,I8)')
' COLVARS| Ring Puckering >>> ring size', &
1332 colvar%ring_puckering_param%nring
1333 IF (colvar%ring_puckering_param%iq == 0)
THEN
1334 WRITE (iw,
'( A,T40,A)')
' COLVARS| Ring Puckering >>> coordinate', &
1335 ' Total Puckering Amplitude'
1336 ELSE IF (colvar%ring_puckering_param%iq > 0)
THEN
1337 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1338 ' Puckering Amplitude', &
1339 colvar%ring_puckering_param%iq
1341 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1342 ' Puckering Angle', &
1343 colvar%ring_puckering_param%iq
1346 WRITE (iw,
'( A)')
' COLVARS| CONDITIONED DISTANCE>> '
1347 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DISTANCE >>> DISTANCE FROM '//tag, &
1348 colvar%mindist_param%i_dist_from(kk), &
1349 kk=1,
SIZE(colvar%mindist_param%i_dist_from))
1350 IF (colvar%mindist_param%use_kinds_from)
THEN
1351 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM KINDS ', &
1352 adjustr(colvar%mindist_param%k_coord_from(kk) (1:10)), &
1353 kk=1,
SIZE(colvar%mindist_param%k_coord_from))
1355 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM '//tag, &
1356 colvar%mindist_param%i_coord_from(kk), &
1357 kk=1,
SIZE(colvar%mindist_param%i_coord_from))
1359 IF (colvar%mindist_param%use_kinds_to)
THEN
1360 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO KINDS ', &
1361 adjustr(colvar%mindist_param%k_coord_to(kk) (1:10)), &
1362 kk=1,
SIZE(colvar%mindist_param%k_coord_to))
1364 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO '//tag, &
1365 colvar%mindist_param%i_coord_to(kk), &
1366 kk=1,
SIZE(colvar%mindist_param%i_coord_to))
1368 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%mindist_param%r_cut
1369 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%mindist_param%p_exp
1370 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%mindist_param%q_exp
1371 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%mindist_param%lambda
1374 IF (colvar%use_points)
THEN
1375 WRITE (iw,
'( A)')
' COLVARS| INFORMATION ON DEFINED GEOMETRICAL POINTS'
1376 DO kk = 1,
SIZE(colvar%points)
1380 WRITE (iw,
'( A)')
' COLVARS| POINT Nr.'//trim(tmpstr2)//
' OF TYPE: '//trim(tmpstr)
1381 IF (
ASSOCIATED(colvar%points(kk)%atoms))
THEN
1382 WRITE (iw,
'( A)')
' COLVARS| ATOMS BUILDING THE GEOMETRICAL POINT'
1383 WRITE (iw,
'( A, I10)') (
' COLVARS| ATOM:', colvar%points(kk)%atoms(k), k=1,
SIZE(colvar%points(kk)%atoms))
1385 WRITE (iw,
'( A,4X,3F12.6)')
' COLVARS| XYZ POSITION OF FIXED POINT:', colvar%points(kk)%r
1391 WRITE (iw,
'( A )')
' '// &
1392 '----------------------------------------------------------------------'
1394 WRITE (iw,
'( A )')
' '// &
1395 '**********************************************************************'
1399 "PRINT%PROGRAM_RUN_INFO")
1400 CALL timestop(handle)
1414 SUBROUTINE read_hydronium_colvars(section, colvar, colvar_id, n_oxygens, n_hydrogens, &
1415 i_oxygens, i_hydrogens)
1418 INTEGER,
INTENT(IN) :: colvar_id
1419 INTEGER,
INTENT(OUT) :: n_oxygens, n_hydrogens
1420 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens, i_hydrogens
1422 INTEGER :: k, n_var, ndim
1423 INTEGER,
DIMENSION(:),
POINTER :: iatms
1431 CALL reallocate(i_oxygens, 1, ndim +
SIZE(iatms))
1432 i_oxygens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1433 ndim = ndim +
SIZE(iatms)
1441 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1442 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1443 ndim = ndim +
SIZE(iatms)
1447 SELECT CASE (colvar_id)
1472 END SUBROUTINE read_hydronium_colvars
1488 SUBROUTINE read_acid_hydronium_colvars(section, colvar, colvar_id, n_oxygens_water, &
1489 n_oxygens_acid, n_hydrogens, i_oxygens_water, &
1490 i_oxygens_acid, i_hydrogens)
1493 INTEGER,
INTENT(IN) :: colvar_id
1494 INTEGER,
INTENT(OUT) :: n_oxygens_water, n_oxygens_acid, &
1496 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens_water, i_oxygens_acid, &
1499 INTEGER :: k, n_var, ndim
1500 INTEGER,
DIMENSION(:),
POINTER :: iatms
1508 CALL reallocate(i_oxygens_water, 1, ndim +
SIZE(iatms))
1509 i_oxygens_water(ndim + 1:ndim +
SIZE(iatms)) = iatms
1510 ndim = ndim +
SIZE(iatms)
1512 n_oxygens_water = ndim
1518 CALL reallocate(i_oxygens_acid, 1, ndim +
SIZE(iatms))
1519 i_oxygens_acid(ndim + 1:ndim +
SIZE(iatms)) = iatms
1520 ndim = ndim +
SIZE(iatms)
1522 n_oxygens_acid = ndim
1528 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1529 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1530 ndim = ndim +
SIZE(iatms)
1534 SELECT CASE (colvar_id)
1565 END SUBROUTINE read_acid_hydronium_colvars
1574 SUBROUTINE colvar_check_points(colvar, section, cell)
1577 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
1579 INTEGER :: i, irep, natoms, npoints, nrep, nweights
1580 INTEGER,
DIMENSION(:),
POINTER :: atoms
1582 REAL(kind=
dp),
DIMENSION(:),
POINTER :: r, weights
1585 NULLIFY (point_sections)
1588 cpassert(
ASSOCIATED(colvar))
1592 colvar%use_points = .true.
1594 ALLOCATE (colvar%points(npoints))
1599 NULLIFY (colvar%points(i)%atoms)
1600 NULLIFY (colvar%points(i)%weights)
1601 CALL section_vals_val_get(point_sections,
"TYPE", i_rep_section=i, i_val=colvar%points(i)%type_id)
1602 SELECT CASE (colvar%points(i)%type_id)
1605 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, n_rep_val=nrep, i_vals=atoms)
1607 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1608 natoms = natoms +
SIZE(atoms)
1610 ALLOCATE (colvar%points(i)%atoms(natoms))
1613 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1614 colvar%points(i)%atoms(natoms + 1:) = atoms(:)
1615 natoms = natoms +
SIZE(atoms)
1618 ALLOCATE (colvar%points(i)%weights(natoms))
1619 colvar%points(i)%weights = 1.0_dp/real(natoms, kind=
dp)
1625 colvar%points(i)%weights(nweights + 1:) = weights(:)
1626 nweights = nweights +
SIZE(weights)
1628 cpassert(natoms == nweights)
1633 colvar%points(i)%r = r
1634 IF (
PRESENT(cell))
THEN
1640 END SUBROUTINE colvar_check_points
1656 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
1659 OPTIONAL,
POINTER :: fixd_list
1662 LOGICAL :: colvar_ok
1664 colvar_ok =
ASSOCIATED(colvar)
1667 IF (
PRESENT(pos))
THEN
1668 DO i = 1,
SIZE(colvar%i_atom)
1669 j = colvar%i_atom(i)
1670 particles(j)%r = pos(:, j)
1674 colvar%dsdr = 0.0_dp
1675 SELECT CASE (colvar%type_id)
1677 CALL dist_colvar(colvar, cell, particles=particles)
1679 CALL coord_colvar(colvar, cell, particles=particles)
1681 CALL population_colvar(colvar, cell, particles=particles)
1683 CALL gyration_radius_colvar(colvar, cell, particles=particles)
1685 CALL torsion_colvar(colvar, cell, particles=particles)
1687 CALL angle_colvar(colvar, cell, particles=particles)
1689 CALL dfunct_colvar(colvar, cell, particles=particles)
1691 CALL plane_distance_colvar(colvar, cell, particles=particles)
1693 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1695 CALL rotation_colvar(colvar, cell, particles=particles)
1697 CALL qparm_colvar(colvar, cell, particles=particles)
1699 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1701 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1703 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1705 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1707 CALL rmsd_colvar(colvar, particles=particles)
1709 CALL reaction_path_colvar(colvar, cell, particles=particles)
1711 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1713 CALL combine_colvar(colvar, cell, particles=particles)
1715 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1717 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1719 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1721 CALL mindist_colvar(colvar, cell, particles=particles)
1723 cpabort(
"need force_env!")
1726 CALL wc_colvar(colvar, cell, particles=particles)
1729 CALL hbp_colvar(colvar, cell, particles=particles)
1731 cpabort(
"Unknown colvar type for colvar_eval_mol_f")
1751 LOGICAL :: colvar_ok
1757 NULLIFY (subsys, cell, colvar, qs_env)
1758 CALL force_env_get(force_env, subsys=subsys, cell=cell, qs_env=qs_env)
1759 colvar_ok =
ASSOCIATED(subsys%colvar_p)
1762 colvar => subsys%colvar_p(icolvar)%colvar
1764 colvar%dsdr = 0.0_dp
1765 SELECT CASE (colvar%type_id)
1767 CALL dist_colvar(colvar, cell, subsys=subsys)
1769 CALL coord_colvar(colvar, cell, subsys=subsys)
1771 CALL population_colvar(colvar, cell, subsys=subsys)
1773 CALL gyration_radius_colvar(colvar, cell, subsys=subsys)
1775 CALL torsion_colvar(colvar, cell, subsys=subsys, no_riemann_sheet_op=.true.)
1777 CALL angle_colvar(colvar, cell, subsys=subsys)
1779 CALL dfunct_colvar(colvar, cell, subsys=subsys)
1781 CALL plane_distance_colvar(colvar, cell, subsys=subsys)
1783 CALL plane_plane_angle_colvar(colvar, cell, subsys=subsys)
1785 CALL rotation_colvar(colvar, cell, subsys=subsys)
1787 CALL qparm_colvar(colvar, cell, subsys=subsys)
1789 CALL hydronium_shell_colvar(colvar, cell, subsys=subsys)
1791 CALL hydronium_dist_colvar(colvar, cell, subsys=subsys)
1793 CALL acid_hyd_dist_colvar(colvar, cell, subsys=subsys)
1795 CALL acid_hyd_shell_colvar(colvar, cell, subsys=subsys)
1797 CALL rmsd_colvar(colvar, subsys=subsys)
1799 CALL reaction_path_colvar(colvar, cell, subsys=subsys)
1801 CALL distance_from_path_colvar(colvar, cell, subsys=subsys)
1803 CALL combine_colvar(colvar, cell, subsys=subsys)
1805 CALL xyz_diag_colvar(colvar, cell, subsys=subsys)
1807 CALL xyz_outerdiag_colvar(colvar, cell, subsys=subsys)
1809 CALL u_colvar(colvar, force_env=force_env)
1811 CALL wc_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1813 CALL hbp_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1815 CALL ring_puckering_colvar(colvar, cell, subsys=subsys)
1817 CALL mindist_colvar(colvar, cell, subsys=subsys)
1819 cpabort(
"Unknown colvar type for colvar_eval_glob_f")
1833 SUBROUTINE colvar_recursive_eval(colvar, cell, particles)
1840 colvar%dsdr = 0.0_dp
1841 SELECT CASE (colvar%type_id)
1843 CALL dist_colvar(colvar, cell, particles=particles)
1845 CALL coord_colvar(colvar, cell, particles=particles)
1847 CALL torsion_colvar(colvar, cell, particles=particles)
1849 CALL angle_colvar(colvar, cell, particles=particles)
1851 CALL dfunct_colvar(colvar, cell, particles=particles)
1853 CALL plane_distance_colvar(colvar, cell, particles=particles)
1855 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1857 CALL rotation_colvar(colvar, cell, particles=particles)
1859 CALL qparm_colvar(colvar, cell, particles=particles)
1861 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1863 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1865 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1867 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1869 CALL rmsd_colvar(colvar, particles=particles)
1871 CALL reaction_path_colvar(colvar, cell, particles=particles)
1873 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1875 CALL combine_colvar(colvar, cell, particles=particles)
1877 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1879 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1881 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1883 CALL mindist_colvar(colvar, cell, particles=particles)
1885 cpabort(
"need force_env!")
1887 CALL wc_colvar(colvar, cell, particles=particles)
1889 CALL hbp_colvar(colvar, cell, particles=particles)
1891 cpabort(
"Unknown colvar type for colvar_recursive_eval")
1893 END SUBROUTINE colvar_recursive_eval
1903 SUBROUTINE get_coordinates(colvar, i, ri, my_particles)
1905 INTEGER,
INTENT(IN) :: i
1906 REAL(kind=
dp),
DIMENSION(3),
INTENT(OUT) :: ri
1909 IF (colvar%use_points)
THEN
1912 ri(:) = my_particles(i)%r(:)
1915 END SUBROUTINE get_coordinates
1925 SUBROUTINE get_mass(colvar, i, mi, my_particles)
1927 INTEGER,
INTENT(IN) :: i
1928 REAL(kind=
dp),
INTENT(OUT) :: mi
1931 IF (colvar%use_points)
THEN
1934 mi = my_particles(i)%atomic_kind%mass
1937 END SUBROUTINE get_mass
1946 SUBROUTINE put_derivative(colvar, i, fi)
1948 INTEGER,
INTENT(IN) :: i
1949 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: fi
1951 IF (colvar%use_points)
THEN
1954 colvar%dsdr(:, i) = colvar%dsdr(:, i) + fi
1957 END SUBROUTINE put_derivative
1967 SUBROUTINE xyz_diag_colvar(colvar, cell, subsys, particles)
1972 POINTER :: particles
1975 REAL(
dp) :: fi(3), r, r0(3), ss(3), xi(3), xpi(3)
1979 NULLIFY (particles_i)
1982 IF (
PRESENT(particles))
THEN
1983 my_particles => particles
1985 cpassert(
PRESENT(subsys))
1987 my_particles => particles_i%els
1989 i = colvar%xyz_diag_param%i_atom
1991 CALL get_coordinates(colvar, i, xpi, my_particles)
1994 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
1995 IF (all(colvar%xyz_diag_param%r0 == huge(0.0_dp)))
THEN
1996 colvar%xyz_diag_param%r0 = xpi
1998 r0 = colvar%xyz_diag_param%r0
2003 IF (colvar%xyz_diag_param%use_pbc)
THEN
2004 ss = matmul(cell%h_inv, xpi - r0)
2006 xi = matmul(cell%hmat, ss)
2011 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
2012 SELECT CASE (colvar%xyz_diag_param%component)
2032 r = xi(1)**2 + xi(2)**2 + xi(3)**2
2035 SELECT CASE (colvar%xyz_diag_param%component)
2052 cpabort(
"xyz_diag_colvar not implemented for anything which is not a single component")
2058 CALL put_derivative(colvar, 1, fi)
2060 END SUBROUTINE xyz_diag_colvar
2070 SUBROUTINE xyz_outerdiag_colvar(colvar, cell, subsys, particles)
2075 POINTER :: particles
2078 REAL(
dp) :: fi(3, 2), r, r0(3), ss(3), xi(3, 2), &
2083 NULLIFY (particles_i)
2086 IF (
PRESENT(particles))
THEN
2087 my_particles => particles
2089 cpassert(
PRESENT(subsys))
2091 my_particles => particles_i%els
2094 i = colvar%xyz_outerdiag_param%i_atoms(k)
2096 CALL get_coordinates(colvar, i, xpi, my_particles)
2097 r0 = colvar%xyz_outerdiag_param%r0(:, k)
2098 IF (all(colvar%xyz_outerdiag_param%r0(:, k) == huge(0.0_dp))) r0 = xpi
2100 IF (colvar%xyz_outerdiag_param%use_pbc)
THEN
2101 ss = matmul(cell%h_inv, xpi - r0)
2103 xi(:, k) = matmul(cell%hmat, ss)
2108 SELECT CASE (colvar%xyz_outerdiag_param%components(k))
2133 IF (xi(l, 1) /= 0.0_dp) fi(l, 1) = fi(l, 1) + xi(i, 2)
2134 r = r + xi(l, 1)*xi(i, 2)
2136 IF (xi(i, 2) /= 0.0_dp) fi(i, 2) = sum(xi(:, 1))
2140 CALL put_derivative(colvar, 1, fi(:, 1))
2141 CALL put_derivative(colvar, 2, fi(:, 2))
2143 END SUBROUTINE xyz_outerdiag_colvar
2153 SUBROUTINE u_colvar(colvar, force_env)
2157 CHARACTER(LEN=default_path_length) :: coupling_function
2158 CHARACTER(LEN=default_string_length) :: def_error, this_error
2159 CHARACTER(LEN=default_string_length), &
2160 DIMENSION(:),
POINTER :: parameters
2161 INTEGER :: iatom, iforce_eval, iparticle, &
2162 jparticle, natom, natom_iforce, &
2164 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms, map_index
2165 REAL(
dp) :: dedf, dx, err, fi(3), lerr, &
2167 REAL(kind=
dp),
DIMENSION(:),
POINTER :: values
2176 IF (
PRESENT(force_env))
THEN
2177 NULLIFY (particles_main, subsys_main)
2179 CALL cp_subsys_get(subsys=subsys_main, particles=particles_main)
2180 natom =
SIZE(particles_main%els)
2181 colvar%n_atom_s = natom
2182 colvar%u_param%natom = natom
2186 colvar%i_atom(iatom) = iatom
2189 IF (.NOT.
ASSOCIATED(colvar%u_param%mixed_energy_section))
THEN
2190 CALL force_env_get(force_env, potential_energy=potential_energy)
2191 colvar%ss = potential_energy
2195 fi(:) = -particles_main%els(iatom)%f
2196 CALL put_derivative(colvar, iatom, fi)
2200 CALL cp_abort(__location__, &
2201 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
2202 ' A combination of mixed force_eval energies has been requested as '// &
2203 ' collective variable, but the MIXED env is not in use! Aborting.')
2205 CALL force_env_get(force_env, force_env_section=force_env_section)
2207 NULLIFY (values, parameters, subsystems, particles, global_forces, map_index, glob_natoms)
2208 nforce_eval =
SIZE(force_env%sub_force_env)
2209 ALLOCATE (glob_natoms(nforce_eval))
2210 ALLOCATE (subsystems(nforce_eval))
2211 ALLOCATE (particles(nforce_eval))
2213 ALLOCATE (global_forces(nforce_eval))
2216 DO iforce_eval = 1, nforce_eval
2217 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
2218 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
2220 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
2221 subsys=subsystems(iforce_eval)%subsys)
2224 particles=particles(iforce_eval)%list)
2227 natom_iforce =
SIZE(particles(iforce_eval)%list%els)
2230 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2231 glob_natoms(iforce_eval) = natom_iforce
2236 CALL force_env%para_env%sync()
2237 CALL force_env%para_env%sum(glob_natoms)
2240 DO iforce_eval = 1, nforce_eval
2241 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
2242 global_forces(iforce_eval)%forces = 0.0_dp
2243 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env))
THEN
2244 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2246 DO iparticle = 1, glob_natoms(iforce_eval)
2247 global_forces(iforce_eval)%forces(:, iparticle) = &
2248 particles(iforce_eval)%list%els(iparticle)%f
2252 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
2255 wrk_section => colvar%u_param%mixed_energy_section
2257 CALL get_generic_info(wrk_section,
"ENERGY_FUNCTION", coupling_function, parameters, &
2258 values, force_env%mixed_env%energies)
2260 CALL parsef(1, trim(coupling_function), parameters)
2262 colvar%ss =
evalf(1, values)
2265 DO iforce_eval = 1, nforce_eval
2268 dedf =
evalfd(1, iforce_eval, values, dx, err)
2269 IF (abs(err) > lerr)
THEN
2270 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
2271 WRITE (def_error,
"(A,G12.6,A)")
"(", lerr,
")"
2274 CALL cp_warn(__location__, &
2275 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
2276 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
2277 trim(def_error)//
' .')
2282 nforce_eval, map_index)
2285 DO iparticle = 1, glob_natoms(iforce_eval)
2286 jparticle = map_index(iparticle)
2287 fi = -dedf*global_forces(iforce_eval)%forces(:, iparticle)
2288 CALL put_derivative(colvar, jparticle, fi)
2291 IF (
ASSOCIATED(map_index))
THEN
2292 DEALLOCATE (map_index)
2296 DO iforce_eval = 1, nforce_eval
2297 DEALLOCATE (global_forces(iforce_eval)%forces)
2299 DEALLOCATE (glob_natoms)
2301 DEALLOCATE (parameters)
2302 DEALLOCATE (global_forces)
2303 DEALLOCATE (subsystems)
2304 DEALLOCATE (particles)
2307 cpabort(
"need force_env!")
2309 END SUBROUTINE u_colvar
2319 SUBROUTINE plane_distance_colvar(colvar, cell, subsys, particles)
2325 POINTER :: particles
2327 INTEGER :: i, j, k, l
2328 REAL(
dp) :: a, b, dsdxpn(3), dxpndxi(3, 3), dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), &
2329 fk(3), fl(3), r12, ri(3), rj(3), rk(3), rl(3), ss(3), xpij(3), xpkj(3), xpl(3), xpn(3)
2333 NULLIFY (particles_i)
2336 IF (
PRESENT(particles))
THEN
2337 my_particles => particles
2339 cpassert(
PRESENT(subsys))
2341 my_particles => particles_i%els
2343 i = colvar%plane_distance_param%plane(1)
2344 j = colvar%plane_distance_param%plane(2)
2345 k = colvar%plane_distance_param%plane(3)
2346 l = colvar%plane_distance_param%point
2348 CALL get_coordinates(colvar, i, ri, my_particles)
2349 CALL get_coordinates(colvar, j, rj, my_particles)
2350 CALL get_coordinates(colvar, k, rk, my_particles)
2351 CALL get_coordinates(colvar, l, rl, my_particles)
2354 xpl = rl - (ri + rj + rk)/3.0_dp
2355 IF (colvar%plane_distance_param%use_pbc)
THEN
2357 ss = matmul(cell%h_inv, ri - rj)
2359 xpij = matmul(cell%hmat, ss)
2361 ss = matmul(cell%h_inv, rk - rj)
2363 xpkj = matmul(cell%hmat, ss)
2365 ss = matmul(cell%h_inv, rl - (ri + rj + rk)/3.0_dp)
2367 xpl = matmul(cell%hmat, ss)
2370 xpn(1) = xpij(2)*xpkj(3) - xpij(3)*xpkj(2)
2371 xpn(2) = xpij(3)*xpkj(1) - xpij(1)*xpkj(3)
2372 xpn(3) = xpij(1)*xpkj(2) - xpij(2)*xpkj(1)
2373 a = dot_product(xpn, xpn)
2374 b = dot_product(xpl, xpn)
2377 dsdxpn(1) = xpl(1)/r12 - b*xpn(1)/(r12*a)
2378 dsdxpn(2) = xpl(2)/r12 - b*xpn(2)/(r12*a)
2379 dsdxpn(3) = xpl(3)/r12 - b*xpn(3)/(r12*a)
2381 dxpndxi(1, 1) = 0.0_dp
2382 dxpndxi(1, 2) = 1.0_dp*xpkj(3)
2383 dxpndxi(1, 3) = -1.0_dp*xpkj(2)
2384 dxpndxi(2, 1) = -1.0_dp*xpkj(3)
2385 dxpndxi(2, 2) = 0.0_dp
2386 dxpndxi(2, 3) = 1.0_dp*xpkj(1)
2387 dxpndxi(3, 1) = 1.0_dp*xpkj(2)
2388 dxpndxi(3, 2) = -1.0_dp*xpkj(1)
2389 dxpndxi(3, 3) = 0.0_dp
2391 dxpndxj(1, 1) = 0.0_dp
2392 dxpndxj(1, 2) = -1.0_dp*xpkj(3) + xpij(3)
2393 dxpndxj(1, 3) = -1.0_dp*xpij(2) + xpkj(2)
2394 dxpndxj(2, 1) = -1.0_dp*xpij(3) + xpkj(3)
2395 dxpndxj(2, 2) = 0.0_dp
2396 dxpndxj(2, 3) = -1.0_dp*xpkj(1) + xpij(1)
2397 dxpndxj(3, 1) = -1.0_dp*xpkj(2) + xpij(2)
2398 dxpndxj(3, 2) = -1.0_dp*xpij(1) + xpkj(1)
2399 dxpndxj(3, 3) = 0.0_dp
2401 dxpndxk(1, 1) = 0.0_dp
2402 dxpndxk(1, 2) = -1.0_dp*xpij(3)
2403 dxpndxk(1, 3) = 1.0_dp*xpij(2)
2404 dxpndxk(2, 1) = 1.0_dp*xpij(3)
2405 dxpndxk(2, 2) = 0.0_dp
2406 dxpndxk(2, 3) = -1.0_dp*xpij(1)
2407 dxpndxk(3, 1) = -1.0_dp*xpij(2)
2408 dxpndxk(3, 2) = 1.0_dp*xpij(1)
2409 dxpndxk(3, 3) = 0.0_dp
2411 fi(:) = matmul(dsdxpn, dxpndxi) - xpn/(3.0_dp*r12)
2412 fj(:) = matmul(dsdxpn, dxpndxj) - xpn/(3.0_dp*r12)
2413 fk(:) = matmul(dsdxpn, dxpndxk) - xpn/(3.0_dp*r12)
2416 CALL put_derivative(colvar, 1, fi)
2417 CALL put_derivative(colvar, 2, fj)
2418 CALL put_derivative(colvar, 3, fk)
2419 CALL put_derivative(colvar, 4, fl)
2421 END SUBROUTINE plane_distance_colvar
2432 SUBROUTINE plane_plane_angle_colvar(colvar, cell, subsys, particles)
2438 POINTER :: particles
2440 INTEGER :: i1, i2, j1, j2, k1, k2, np
2442 REAL(
dp) :: a1, a2, d, dnorm_dxpn(3), dprod12_dxpn(3), dsdxpn(3), dt_dxpn(3), dxpndxi(3, 3), &
2443 dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), fk(3), fmod, norm1, norm2, prod_12, ri1(3), &
2444 ri2(3), rj1(3), rj2(3), rk1(3), rk2(3), ss(3), t, xpij1(3), xpij2(3), xpkj1(3), xpkj2(3), &
2449 NULLIFY (particles_i)
2453 IF (
PRESENT(particles))
THEN
2454 my_particles => particles
2456 cpassert(
PRESENT(subsys))
2458 my_particles => particles_i%els
2462 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2463 i1 = colvar%plane_plane_angle_param%plane1%points(1)
2464 j1 = colvar%plane_plane_angle_param%plane1%points(2)
2465 k1 = colvar%plane_plane_angle_param%plane1%points(3)
2468 CALL get_coordinates(colvar, i1, ri1, my_particles)
2469 CALL get_coordinates(colvar, j1, rj1, my_particles)
2470 CALL get_coordinates(colvar, k1, rk1, my_particles)
2473 ss = matmul(cell%h_inv, ri1 - rj1)
2475 xpij1 = matmul(cell%hmat, ss)
2478 ss = matmul(cell%h_inv, rk1 - rj1)
2480 xpkj1 = matmul(cell%hmat, ss)
2483 xpn1(1) = xpij1(2)*xpkj1(3) - xpij1(3)*xpkj1(2)
2484 xpn1(2) = xpij1(3)*xpkj1(1) - xpij1(1)*xpkj1(3)
2485 xpn1(3) = xpij1(1)*xpkj1(2) - xpij1(2)*xpkj1(1)
2487 xpn1 = colvar%plane_plane_angle_param%plane1%normal_vec
2489 a1 = dot_product(xpn1, xpn1)
2491 cpassert(norm1 /= 0.0_dp)
2494 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2495 i2 = colvar%plane_plane_angle_param%plane2%points(1)
2496 j2 = colvar%plane_plane_angle_param%plane2%points(2)
2497 k2 = colvar%plane_plane_angle_param%plane2%points(3)
2500 CALL get_coordinates(colvar, i2, ri2, my_particles)
2501 CALL get_coordinates(colvar, j2, rj2, my_particles)
2502 CALL get_coordinates(colvar, k2, rk2, my_particles)
2505 ss = matmul(cell%h_inv, ri2 - rj2)
2507 xpij2 = matmul(cell%hmat, ss)
2510 ss = matmul(cell%h_inv, rk2 - rj2)
2512 xpkj2 = matmul(cell%hmat, ss)
2515 xpn2(1) = xpij2(2)*xpkj2(3) - xpij2(3)*xpkj2(2)
2516 xpn2(2) = xpij2(3)*xpkj2(1) - xpij2(1)*xpkj2(3)
2517 xpn2(3) = xpij2(1)*xpkj2(2) - xpij2(2)*xpkj2(1)
2519 xpn2 = colvar%plane_plane_angle_param%plane2%normal_vec
2521 a2 = dot_product(xpn2, xpn2)
2523 cpassert(norm2 /= 0.0_dp)
2526 prod_12 = dot_product(xpn1, xpn2)
2530 t = min(1.0_dp, abs(t))*sign(1.0_dp, t)
2533 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2536 fmod = -1.0_dp/sin(colvar%ss)
2541 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2543 dnorm_dxpn = 1.0_dp/norm1*xpn1
2544 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm2)/d**2
2546 dsdxpn(1) = fmod*dt_dxpn(1)
2547 dsdxpn(2) = fmod*dt_dxpn(2)
2548 dsdxpn(3) = fmod*dt_dxpn(3)
2550 dxpndxi(1, 1) = 0.0_dp
2551 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2552 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2553 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2554 dxpndxi(2, 2) = 0.0_dp
2555 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2556 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2557 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2558 dxpndxi(3, 3) = 0.0_dp
2560 dxpndxj(1, 1) = 0.0_dp
2561 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2562 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2563 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2564 dxpndxj(2, 2) = 0.0_dp
2565 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2566 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2567 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2568 dxpndxj(3, 3) = 0.0_dp
2570 dxpndxk(1, 1) = 0.0_dp
2571 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2572 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2573 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2574 dxpndxk(2, 2) = 0.0_dp
2575 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2576 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2577 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2578 dxpndxk(3, 3) = 0.0_dp
2580 fi = matmul(dsdxpn, dxpndxi)
2581 fj = matmul(dsdxpn, dxpndxj)
2582 fk = matmul(dsdxpn, dxpndxk)
2585 CALL put_derivative(colvar, np + 1, fi)
2586 CALL put_derivative(colvar, np + 2, fj)
2587 CALL put_derivative(colvar, np + 3, fk)
2592 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2594 dnorm_dxpn = 1.0_dp/norm2*xpn2
2595 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm1)/d**2
2597 dsdxpn(1) = fmod*dt_dxpn(1)
2598 dsdxpn(2) = fmod*dt_dxpn(2)
2599 dsdxpn(3) = fmod*dt_dxpn(3)
2601 dxpndxi(1, 1) = 0.0_dp
2602 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2603 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2604 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2605 dxpndxi(2, 2) = 0.0_dp
2606 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2607 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2608 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2609 dxpndxi(3, 3) = 0.0_dp
2611 dxpndxj(1, 1) = 0.0_dp
2612 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2613 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2614 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2615 dxpndxj(2, 2) = 0.0_dp
2616 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2617 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2618 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2619 dxpndxj(3, 3) = 0.0_dp
2621 dxpndxk(1, 1) = 0.0_dp
2622 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2623 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2624 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2625 dxpndxk(2, 2) = 0.0_dp
2626 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2627 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2628 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2629 dxpndxk(3, 3) = 0.0_dp
2631 fi = matmul(dsdxpn, dxpndxi)
2632 fj = matmul(dsdxpn, dxpndxj)
2633 fk = matmul(dsdxpn, dxpndxk)
2636 CALL put_derivative(colvar, np + 1, fi)
2637 CALL put_derivative(colvar, np + 2, fj)
2638 CALL put_derivative(colvar, np + 3, fk)
2641 END SUBROUTINE plane_plane_angle_colvar
2651 SUBROUTINE rotation_colvar(colvar, cell, subsys, particles)
2656 POINTER :: particles
2659 REAL(
dp) :: a, b, fmod, t0, t1, t2, t3, xdum(3), &
2661 REAL(kind=
dp) :: dp1b1(3), dp1b2(3), dp2b1(3), dp2b2(3), &
2662 ss(3), xp1b1(3), xp1b2(3), xp2b1(3), &
2667 NULLIFY (particles_i)
2670 IF (
PRESENT(particles))
THEN
2671 my_particles => particles
2673 cpassert(
PRESENT(subsys))
2675 my_particles => particles_i%els
2677 i = colvar%rotation_param%i_at1_bond1
2678 CALL get_coordinates(colvar, i, xp1b1, my_particles)
2679 i = colvar%rotation_param%i_at2_bond1
2680 CALL get_coordinates(colvar, i, xp2b1, my_particles)
2681 i = colvar%rotation_param%i_at1_bond2
2682 CALL get_coordinates(colvar, i, xp1b2, my_particles)
2683 i = colvar%rotation_param%i_at2_bond2
2684 CALL get_coordinates(colvar, i, xp2b2, my_particles)
2686 ss = matmul(cell%h_inv, xp1b1 - xp2b1)
2688 xij = matmul(cell%hmat, ss)
2690 ss = matmul(cell%h_inv, xp1b2 - xp2b2)
2692 xkj = matmul(cell%hmat, ss)
2697 t1 = 1.0_dp/(a**3.0_dp*b)
2698 t2 = 1.0_dp/(a*b**3.0_dp)
2699 t3 = dot_product(xij, xkj)
2700 colvar%ss = acos(t3*t0)
2701 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2704 fmod = -1.0_dp/sin(colvar%ss)
2706 dp1b1 = xkj(:)*t0 - xij(:)*t1*t3
2707 dp2b1 = -xkj(:)*t0 + xij(:)*t1*t3
2708 dp1b2 = xij(:)*t0 - xkj(:)*t2*t3
2709 dp2b2 = -xij(:)*t0 + xkj(:)*t2*t3
2712 idum = colvar%rotation_param%i_at1_bond1
2713 CALL put_derivative(colvar, idum, xdum)
2715 idum = colvar%rotation_param%i_at2_bond1
2716 CALL put_derivative(colvar, idum, xdum)
2718 idum = colvar%rotation_param%i_at1_bond2
2719 CALL put_derivative(colvar, idum, xdum)
2721 idum = colvar%rotation_param%i_at2_bond2
2722 CALL put_derivative(colvar, idum, xdum)
2724 END SUBROUTINE rotation_colvar
2735 SUBROUTINE dfunct_colvar(colvar, cell, subsys, particles)
2740 POINTER :: particles
2742 INTEGER :: i, j, k, l
2743 REAL(
dp) :: fi(3), fj(3), fk(3), fl(3), r12, r34, &
2744 ss(3), xij(3), xkl(3), xpi(3), xpj(3), &
2749 NULLIFY (particles_i)
2752 IF (
PRESENT(particles))
THEN
2753 my_particles => particles
2755 cpassert(
PRESENT(subsys))
2757 my_particles => particles_i%els
2759 i = colvar%dfunct_param%i_at_dfunct(1)
2760 j = colvar%dfunct_param%i_at_dfunct(2)
2762 CALL get_coordinates(colvar, i, xpi, my_particles)
2763 CALL get_coordinates(colvar, j, xpj, my_particles)
2764 IF (colvar%dfunct_param%use_pbc)
THEN
2765 ss = matmul(cell%h_inv, xpi - xpj)
2767 xij = matmul(cell%hmat, ss)
2771 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2773 k = colvar%dfunct_param%i_at_dfunct(3)
2774 l = colvar%dfunct_param%i_at_dfunct(4)
2775 CALL get_coordinates(colvar, k, xpk, my_particles)
2776 CALL get_coordinates(colvar, l, xpl, my_particles)
2777 IF (colvar%dfunct_param%use_pbc)
THEN
2778 ss = matmul(cell%h_inv, xpk - xpl)
2780 xkl = matmul(cell%hmat, ss)
2784 r34 = sqrt(xkl(1)**2 + xkl(2)**2 + xkl(3)**2)
2786 colvar%ss = r12 + colvar%dfunct_param%coeff*r34
2789 fk(:) = colvar%dfunct_param%coeff*xkl/r34
2790 fl(:) = -colvar%dfunct_param%coeff*xkl/r34
2791 CALL put_derivative(colvar, 1, fi)
2792 CALL put_derivative(colvar, 2, fj)
2793 CALL put_derivative(colvar, 3, fk)
2794 CALL put_derivative(colvar, 4, fl)
2796 END SUBROUTINE dfunct_colvar
2806 SUBROUTINE angle_colvar(colvar, cell, subsys, particles)
2811 POINTER :: particles
2814 REAL(
dp) :: a, b, fi(3), fj(3), fk(3), fmod, ri(3), &
2815 rj(3), rk(3), ss(3), t0, t1, t2, t3, &
2820 NULLIFY (particles_i)
2823 IF (
PRESENT(particles))
THEN
2824 my_particles => particles
2826 cpassert(
PRESENT(subsys))
2828 my_particles => particles_i%els
2830 i = colvar%angle_param%i_at_angle(1)
2831 j = colvar%angle_param%i_at_angle(2)
2832 k = colvar%angle_param%i_at_angle(3)
2833 CALL get_coordinates(colvar, i, ri, my_particles)
2834 CALL get_coordinates(colvar, j, rj, my_particles)
2835 CALL get_coordinates(colvar, k, rk, my_particles)
2837 ss = matmul(cell%h_inv, ri - rj)
2839 xij = matmul(cell%hmat, ss)
2841 ss = matmul(cell%h_inv, rk - rj)
2843 xkj = matmul(cell%hmat, ss)
2848 t1 = 1.0_dp/(a**3.0_dp*b)
2849 t2 = 1.0_dp/(a*b**3.0_dp)
2850 t3 = dot_product(xij, xkj)
2851 colvar%ss = acos(t3*t0)
2852 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2855 fmod = -1.0_dp/sin(colvar%ss)
2857 fi(:) = xkj(:)*t0 - xij(:)*t1*t3
2858 fj(:) = -xkj(:)*t0 + xij(:)*t1*t3 - xij(:)*t0 + xkj(:)*t2*t3
2859 fk(:) = xij(:)*t0 - xkj(:)*t2*t3
2863 CALL put_derivative(colvar, 1, fi)
2864 CALL put_derivative(colvar, 2, fj)
2865 CALL put_derivative(colvar, 3, fk)
2867 END SUBROUTINE angle_colvar
2877 SUBROUTINE dist_colvar(colvar, cell, subsys, particles)
2882 POINTER :: particles
2885 REAL(
dp) :: fi(3), fj(3), r12, ss(3), xij(3), &
2890 NULLIFY (particles_i)
2893 IF (
PRESENT(particles))
THEN
2894 my_particles => particles
2896 cpassert(
PRESENT(subsys))
2898 my_particles => particles_i%els
2900 i = colvar%dist_param%i_at
2901 j = colvar%dist_param%j_at
2902 CALL get_coordinates(colvar, i, xpi, my_particles)
2903 CALL get_coordinates(colvar, j, xpj, my_particles)
2904 ss = matmul(cell%h_inv, xpi - xpj)
2906 xij = matmul(cell%hmat, ss)
2907 SELECT CASE (colvar%dist_param%axis_id)
2926 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2928 IF (colvar%dist_param%sign_d)
THEN
2929 SELECT CASE (colvar%dist_param%axis_id)
2947 CALL put_derivative(colvar, 1, fi)
2948 CALL put_derivative(colvar, 2, fj)
2950 END SUBROUTINE dist_colvar
2961 SUBROUTINE torsion_colvar(colvar, cell, subsys, particles, no_riemann_sheet_op)
2967 POINTER :: particles
2968 LOGICAL,
INTENT(IN),
OPTIONAL :: no_riemann_sheet_op
2971 LOGICAL :: no_riemann_sheet
2972 REAL(
dp) ::
angle, cosine, dedphi, dedxia, dedxib, dedxic, dedxid, dedxt, dedxu, dedyia, &
2973 dedyib, dedyic, dedyid, dedyt, dedyu, dedzia, dedzib, dedzic, dedzid, dedzt, dedzu, dt, &
2974 e, ftmp(3), o0, rcb, rt2, rtmp(3), rtru, ru2, sine, ss(3), xba, xca, xcb, xdb, xdc, xt, &
2975 xtu, xu, yba, yca, ycb, ydb, ydc, yt, ytu, yu, zba, zca, zcb, zdb, zdc, zt, ztu, zu
2976 REAL(
dp),
DIMENSION(3, 4) :: rr
2980 NULLIFY (particles_i)
2982 IF (
PRESENT(particles))
THEN
2983 my_particles => particles
2985 cpassert(
PRESENT(subsys))
2987 my_particles => particles_i%els
2989 no_riemann_sheet = .false.
2990 IF (
PRESENT(no_riemann_sheet_op)) no_riemann_sheet = no_riemann_sheet_op
2992 i = colvar%torsion_param%i_at_tors(ii)
2993 CALL get_coordinates(colvar, i, rtmp, my_particles)
2994 rr(:, ii) = rtmp(1:3)
2996 o0 = colvar%torsion_param%o0
2998 ss = matmul(cell%h_inv, rr(:, 2) - rr(:, 1))
3000 ss = matmul(cell%hmat, ss)
3005 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 2))
3007 ss = matmul(cell%hmat, ss)
3012 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 3))
3014 ss = matmul(cell%hmat, ss)
3019 xt = yba*zcb - ycb*zba
3020 yt = zba*xcb - zcb*xba
3021 zt = xba*ycb - xcb*yba
3022 xu = ycb*zdc - ydc*zcb
3023 yu = zcb*xdc - zdc*xcb
3024 zu = xcb*ydc - xdc*ycb
3028 rt2 = xt*xt + yt*yt + zt*zt
3029 ru2 = xu*xu + yu*yu + zu*zu
3030 rtru = sqrt(rt2*ru2)
3031 IF (rtru /= 0.0_dp)
THEN
3032 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb)
3033 cosine = (xt*xu + yt*yu + zt*zu)/rtru
3034 sine = (xcb*xtu + ycb*ytu + zcb*ztu)/(rcb*rtru)
3035 cosine = min(1.0_dp, max(-1.0_dp, cosine))
3036 angle = acos(cosine)
3040 dt = mod(2.0e4_dp*
pi + dt - o0, 2.0_dp*
pi)
3041 IF (dt >
pi) dt = dt - 2.0_dp*
pi
3043 colvar%torsion_param%o0 = dt
3053 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 1))
3055 ss = matmul(cell%hmat, ss)
3060 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 2))
3062 ss = matmul(cell%hmat, ss)
3067 dedxt = dedphi*(yt*zcb - ycb*zt)/(rt2*rcb)
3068 dedyt = dedphi*(zt*xcb - zcb*xt)/(rt2*rcb)
3069 dedzt = dedphi*(xt*ycb - xcb*yt)/(rt2*rcb)
3070 dedxu = -dedphi*(yu*zcb - ycb*zu)/(ru2*rcb)
3071 dedyu = -dedphi*(zu*xcb - zcb*xu)/(ru2*rcb)
3072 dedzu = -dedphi*(xu*ycb - xcb*yu)/(ru2*rcb)
3076 dedxia = zcb*dedyt - ycb*dedzt
3077 dedyia = xcb*dedzt - zcb*dedxt
3078 dedzia = ycb*dedxt - xcb*dedyt
3079 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu
3080 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu
3081 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu
3082 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu
3083 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu
3084 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu
3085 dedxid = zcb*dedyu - ycb*dedzu
3086 dedyid = xcb*dedzu - zcb*dedxu
3087 dedzid = ycb*dedxu - xcb*dedyu
3104 IF (no_riemann_sheet) colvar%ss = atan2(sin(e), cos(e))
3108 CALL put_derivative(colvar, 1, ftmp)
3112 CALL put_derivative(colvar, 2, ftmp)
3116 CALL put_derivative(colvar, 3, ftmp)
3120 CALL put_derivative(colvar, 4, ftmp)
3121 END SUBROUTINE torsion_colvar
3130 SUBROUTINE qparm_colvar(colvar, cell, subsys, particles)
3135 POINTER :: particles
3137 INTEGER :: aa, bb, cc, i, idim, ii, j, jj, l, mm, &
3138 n_atoms_from, n_atoms_to, ncells(3)
3139 LOGICAL :: include_images
3140 REAL(kind=
dp) :: denominator_tolerance, fact, ftmp(3), im_qlm, inv_n_atoms_from, nbond, &
3141 pre_fac, ql, qparm, r1cut, rcut, re_qlm, rij, rij_shift, shift(3), ss(3), ss0(3), xij(3), &
3143 REAL(kind=
dp),
DIMENSION(3) :: d_im_qlm_dxi, d_nbond_dxi, d_ql_dxi, &
3144 d_re_qlm_dxi, xpi, xpj
3148 n_atoms_to = colvar%qparm_param%n_atoms_to
3149 n_atoms_from = colvar%qparm_param%n_atoms_from
3150 rcut = colvar%qparm_param%rcut
3151 l = colvar%qparm_param%l
3152 r1cut = colvar%qparm_param%rstart
3153 include_images = colvar%qparm_param%include_images
3154 NULLIFY (particles_i)
3156 IF (
PRESENT(particles))
THEN
3157 my_particles => particles
3159 cpassert(
PRESENT(subsys))
3161 my_particles => particles_i%els
3163 cpassert(r1cut < rcut)
3164 denominator_tolerance = 1.0e-8_dp
3167 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
3168 DO ii = 1, n_atoms_from
3169 i = colvar%qparm_param%i_at_from(ii)
3170 CALL get_coordinates(colvar, i, xpi, my_particles)
3173 d_ql_dxi(:) = 0.0_dp
3179 d_re_qlm_dxi(:) = 0.0_dp
3180 d_im_qlm_dxi(:) = 0.0_dp
3181 d_nbond_dxi(:) = 0.0_dp
3183 jloop:
DO jj = 1, n_atoms_to
3185 j = colvar%qparm_param%i_at_to(jj)
3186 CALL get_coordinates(colvar, j, xpj, my_particles)
3188 IF (include_images)
THEN
3190 cpassert(cell%orthorhombic)
3194 xij(:) = xpj(:) - xpi(:)
3195 ss = matmul(cell%h_inv, xij)
3201 shift(idim) = 1.0_dp
3202 xij_shift = matmul(cell%hmat, shift)
3203 rij_shift = norm2(xij_shift)
3204 ncells(idim) = floor(rcut/rij_shift - 0.5)
3208 DO aa = -ncells(1), ncells(1)
3209 DO bb = -ncells(2), ncells(2)
3210 DO cc = -ncells(3), ncells(3)
3212 IF (i == j .AND. aa == 0 .AND. bb == 0 .AND. cc == 0) cycle
3213 shift(1) = real(aa, kind=
dp)
3214 shift(2) = real(bb, kind=
dp)
3215 shift(3) = real(cc, kind=
dp)
3216 xij = matmul(cell%hmat, ss0(:) + shift(:))
3218 IF (rij > rcut) cycle
3221 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3222 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3223 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3231 IF (i == j) cycle jloop
3232 xij(:) = xpj(:) - xpi(:)
3234 IF (rij > rcut) cycle jloop
3237 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3238 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3239 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3249 IF (nbond < denominator_tolerance)
THEN
3250 cpwarn(
"QPARM: number of neighbors is very close to zero!")
3253 d_nbond_dxi(:) = d_nbond_dxi(:)/nbond
3254 re_qlm = re_qlm/nbond
3255 d_re_qlm_dxi(:) = d_re_qlm_dxi(:)/nbond - d_nbond_dxi(:)*re_qlm
3256 im_qlm = im_qlm/nbond
3257 d_im_qlm_dxi(:) = d_im_qlm_dxi(:)/nbond - d_nbond_dxi(:)*im_qlm
3259 ql = ql + fact*(re_qlm*re_qlm + im_qlm*im_qlm)
3260 d_ql_dxi(:) = d_ql_dxi(:) &
3261 + fact*2.0_dp*(re_qlm*d_re_qlm_dxi(:) + im_qlm*d_im_qlm_dxi(:))
3265 pre_fac = (4.0_dp*
pi)/(2.0_dp*l + 1)
3266 qparm = qparm + sqrt(pre_fac*ql)
3267 ftmp(:) = 0.5_dp*sqrt(pre_fac/ql)*d_ql_dxi(:)
3269 ftmp(:) = -1.0_dp*ftmp(:)
3271 CALL put_derivative(colvar, ii, ftmp)
3275 colvar%ss = qparm*inv_n_atoms_from
3276 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
3278 END SUBROUTINE qparm_colvar
3296 SUBROUTINE accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3297 denominator_tolerance, ll, mm, nbond, re_qlm, im_qlm, &
3298 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3300 REAL(kind=
dp),
INTENT(IN) :: xij(3), rij, rcut, r1cut, &
3301 denominator_tolerance
3302 INTEGER,
INTENT(IN) :: ll, mm
3303 REAL(kind=
dp),
INTENT(INOUT) :: nbond, re_qlm, im_qlm, d_re_qlm_dxi(3), &
3304 d_im_qlm_dxi(3), d_nbond_dxi(3)
3306 REAL(kind=
dp) :: bond, costheta, dplm, dylm, exp0, &
3307 exp_fac, fi, plm, pre_fac, sqrt_c1
3308 REAL(kind=
dp),
DIMENSION(3) :: dcostheta, dfi
3312 IF (rij > rcut)
THEN
3317 IF (rij < r1cut)
THEN
3321 exp0 = exp((r1cut - rcut)/(rij - rcut) - (r1cut - rcut)/(r1cut - rij))
3322 bond = 1.0_dp/(1.0_dp + exp0)
3323 exp_fac = ((rcut - r1cut)/(rij - rcut)**2 + (rcut - r1cut)/(r1cut - rij)**2)*exp0/(1.0_dp + exp0)**2
3326 IF (bond > 1.0_dp)
THEN
3327 cpabort(
"bond > 1.0_dp")
3330 nbond = nbond + bond
3331 IF (abs(xij(1)) < denominator_tolerance &
3332 .AND. abs(xij(2)) < denominator_tolerance)
THEN
3335 fi = atan2(xij(2), xij(1))
3338 costheta = xij(3)/rij
3339 IF (costheta > 1.0_dp) costheta = 1.0_dp
3340 IF (costheta < -1.0_dp) costheta = -1.0_dp
3345 IF ((ll + abs(mm)) >
maxfac)
THEN
3346 cpabort(
"(l+m) > maxfac")
3349 sqrt_c1 = sqrt(((2*ll + 1)*
fac(ll - abs(mm)))/(4*
pi*
fac(ll + abs(mm))))
3350 pre_fac = bond*sqrt_c1
3353 re_qlm = re_qlm + pre_fac*plm*cos(mm*fi)
3354 im_qlm = im_qlm + pre_fac*plm*sin(mm*fi)
3356 dcostheta(:) = xij(:)*xij(3)/(rij**3)
3357 dcostheta(3) = dcostheta(3) - 1.0_dp/rij
3361 dfi(1) = xij(2)/(xij(1)**2 + xij(2)**2)
3362 dfi(2) = -xij(1)/(xij(1)**2 + xij(2)**2)
3364 d_re_qlm_dxi(:) = d_re_qlm_dxi(:) &
3365 + exp_fac*sqrt_c1*plm*cos(mm*fi)*xij(:)/rij &
3366 + dylm*dcostheta(:)*cos(mm*fi) &
3367 + pre_fac*plm*mm*(-1.0_dp)*sin(mm*fi)*dfi(:)
3368 d_im_qlm_dxi(:) = d_im_qlm_dxi(:) &
3369 + exp_fac*sqrt_c1*plm*sin(mm*fi)*xij(:)/rij &
3370 + dylm*dcostheta(:)*sin(mm*fi) &
3371 + pre_fac*plm*mm*(+1.0_dp)*cos(mm*fi)*dfi(:)
3372 d_nbond_dxi(:) = d_nbond_dxi(:) + exp_fac*xij(:)/rij
3374 END SUBROUTINE accumulate_qlm_over_neigbors
3386 SUBROUTINE hydronium_shell_colvar(colvar, cell, subsys, particles)
3391 POINTER :: particles
3393 INTEGER :: i, ii, j, jj, n_hydrogens, n_oxygens, &
3394 pm, poh, poo, qm, qoh, qoo
3395 REAL(
dp) :: drji, fscalar, invden, lambda, nh, num, &
3396 qtot, rji(3), roh, roo, rrel
3397 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: m, noh, noo, qloc
3398 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dm, dnoh, dnoo
3399 REAL(
dp),
DIMENSION(3) :: rpi, rpj
3403 n_oxygens = colvar%hydronium_shell_param%n_oxygens
3404 n_hydrogens = colvar%hydronium_shell_param%n_hydrogens
3405 nh = colvar%hydronium_shell_param%nh
3406 poh = colvar%hydronium_shell_param%poh
3407 qoh = colvar%hydronium_shell_param%qoh
3408 poo = colvar%hydronium_shell_param%poo
3409 qoo = colvar%hydronium_shell_param%qoo
3410 roo = colvar%hydronium_shell_param%roo
3411 roh = colvar%hydronium_shell_param%roh
3412 lambda = colvar%hydronium_shell_param%lambda
3413 pm = colvar%hydronium_shell_param%pm
3414 qm = colvar%hydronium_shell_param%qm
3416 NULLIFY (particles_i)
3418 IF (
PRESENT(particles))
THEN
3419 my_particles => particles
3421 cpassert(
PRESENT(subsys))
3423 my_particles => particles_i%els
3426 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3427 ALLOCATE (noh(n_oxygens))
3428 ALLOCATE (m(n_oxygens))
3429 ALLOCATE (dm(3, n_hydrogens, n_oxygens))
3431 ALLOCATE (dnoo(3, n_oxygens, n_oxygens))
3432 ALLOCATE (noo(n_oxygens))
3434 ALLOCATE (qloc(n_oxygens))
3444 DO ii = 1, n_oxygens
3445 i = colvar%hydronium_shell_param%i_oxygens(ii)
3446 rpi(:) = my_particles(i)%r(1:3)
3448 DO jj = 1, n_hydrogens
3449 j = colvar%hydronium_shell_param%i_hydrogens(jj)
3450 rpj(:) = my_particles(j)%r(1:3)
3451 rji =
pbc(rpj, rpi, cell)
3452 drji = sqrt(sum(rji**2))
3454 num = (1.0_dp - rrel**poh)
3455 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3456 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3457 noh(ii) = noh(ii) + num*invden
3458 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3459 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3460 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3463 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3464 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3465 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3468 m(ii) = 1.0_dp - (1.0_dp - (noh(ii)/nh)**pm)/ &
3469 (1.0_dp - (noh(ii)/nh)**qm)
3472 DO jj = 1, n_oxygens
3474 j = colvar%hydronium_shell_param%i_oxygens(jj)
3475 rpj(:) = my_particles(j)%r(1:3)
3476 rji =
pbc(rpj, rpi, cell)
3477 drji = sqrt(sum(rji**2))
3479 num = (1.0_dp - rrel**poo)
3480 invden = 1.0_dp/(1.0_dp - rrel**qoo)
3481 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3482 noo(ii) = noo(ii) + num*invden
3483 fscalar = ((-poo*(rrel**(poo - 1))*invden) &
3484 + num*(invden)**2*qoo*(rrel**(qoo - 1)))/(drji*roo)
3485 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3488 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
3489 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drji)
3490 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3497 DO ii = 1, n_oxygens
3498 qloc(ii) = exp(lambda*m(ii)*noo(ii))
3499 qtot = qtot + qloc(ii)
3502 DO ii = 1, n_oxygens
3504 DO jj = 1, n_hydrogens
3505 dm(1:3, jj, ii) = (pm*((noh(ii)/nh)**(pm - 1))*dnoh(1:3, jj, ii))/nh/ &
3506 (1.0_dp - (noh(ii)/nh)**qm) - &
3507 (1.0_dp - (noh(ii)/nh)**pm)/ &
3508 ((1.0_dp - (noh(ii)/nh)**qm)**2)* &
3509 qm*dnoh(1:3, jj, ii)*(noh(ii)/nh)**(qm - 1)/nh
3511 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3512 colvar%dsdr(1:3, n_oxygens + jj) = colvar%dsdr(1:3, n_oxygens + jj) &
3513 - qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3516 DO jj = 1, n_oxygens
3517 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3518 colvar%dsdr(1:3, jj) = colvar%dsdr(1:3, jj) &
3519 - qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3523 colvar%ss = log(qtot)/lambda
3532 END SUBROUTINE hydronium_shell_colvar
3543 SUBROUTINE hydronium_dist_colvar(colvar, cell, subsys, particles)
3548 POINTER :: particles
3550 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3551 n_oxygens, offseth, pf, pm, poh, qf, &
3553 REAL(
dp) :: drji, drki, fscalar, invden, lambda, nh, nn, num, rion, rion_den, rion_num, &
3554 rji(3), rki(3), roh, rrel, sum_expfac_f, sum_expfac_noh
3555 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac_f, dexpfac_noh, df, dm, &
3556 expfac_f, expfac_f_rki, expfac_noh, f, &
3558 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_f_rki
3559 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rki, dnoh
3560 REAL(
dp),
DIMENSION(3) :: rpi, rpj, rpk
3564 n_oxygens = colvar%hydronium_dist_param%n_oxygens
3565 n_hydrogens = colvar%hydronium_dist_param%n_hydrogens
3566 poh = colvar%hydronium_dist_param%poh
3567 qoh = colvar%hydronium_dist_param%qoh
3568 roh = colvar%hydronium_dist_param%roh
3569 pm = colvar%hydronium_dist_param%pm
3570 qm = colvar%hydronium_dist_param%qm
3571 nh = colvar%hydronium_dist_param%nh
3572 pf = colvar%hydronium_dist_param%pf
3573 qf = colvar%hydronium_dist_param%qf
3574 nn = colvar%hydronium_dist_param%nn
3575 lambda = colvar%hydronium_dist_param%lambda
3577 NULLIFY (particles_i)
3579 IF (
PRESENT(particles))
THEN
3580 my_particles => particles
3582 cpassert(
PRESENT(subsys))
3584 my_particles => particles_i%els
3587 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3588 ALLOCATE (noh(n_oxygens))
3589 ALLOCATE (m(n_oxygens), dm(n_oxygens))
3590 ALLOCATE (f(n_oxygens), df(n_oxygens))
3591 ALLOCATE (expfac_noh(n_oxygens), dexpfac_noh(n_oxygens))
3592 ALLOCATE (expfac_f(n_oxygens), dexpfac_f(n_oxygens))
3593 ALLOCATE (ddist_rki(3, n_oxygens, n_oxygens))
3594 ALLOCATE (expfac_f_rki(n_oxygens))
3595 ALLOCATE (dexpfac_f_rki(n_oxygens, n_oxygens))
3607 sum_expfac_noh = 0._dp
3608 sum_expfac_f = 0._dp
3610 expfac_f_rki = 0._dp
3611 dexpfac_f_rki = 0._dp
3614 DO ii = 1, n_oxygens
3615 i = colvar%hydronium_dist_param%i_oxygens(ii)
3616 rpi(:) = my_particles(i)%r(1:3)
3617 DO jj = 1, n_hydrogens
3618 j = colvar%hydronium_dist_param%i_hydrogens(jj)
3619 rpj(:) = my_particles(j)%r(1:3)
3620 rji =
pbc(rpj, rpi, cell)
3621 drji = sqrt(sum(rji**2))
3623 num = (1.0_dp - rrel**poh)
3624 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3625 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3626 noh(ii) = noh(ii) + num*invden
3627 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3628 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3629 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3632 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3633 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3634 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3640 DO ii = 1, n_oxygens
3641 num = 1.0_dp - (noh(ii)/nh)**pm
3642 invden = 1.0_dp/(1.0_dp - (noh(ii)/nh)**qm)
3643 m(ii) = 1.0_dp - num*invden
3644 dm(ii) = (pm*(noh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
3645 (noh(ii)/nh)**(qm - 1))/nh
3646 expfac_noh(ii) = exp(lambda*noh(ii))
3647 dexpfac_noh(ii) = lambda*expfac_noh(ii)
3648 sum_expfac_noh = sum_expfac_noh + expfac_noh(ii)
3652 DO ii = 1, n_oxygens
3653 i = colvar%hydronium_dist_param%i_oxygens(ii)
3654 num = 1.0_dp - (noh(ii)/nn)**pf
3655 invden = 1.0_dp/(1.0_dp - (noh(ii)/nn)**qf)
3657 df(ii) = (-pf*(noh(ii)/nn)**(pf - 1)*invden + qf*num*(invden**2)* &
3658 (noh(ii)/nn)**(qf - 1))/nn
3659 expfac_f(ii) = exp(lambda*f(ii))
3660 dexpfac_f(ii) = lambda*expfac_f(ii)
3661 sum_expfac_f = sum_expfac_f + expfac_f(ii)
3665 DO ii = 1, n_oxygens
3666 i = colvar%hydronium_dist_param%i_oxygens(ii)
3667 rpi(:) = my_particles(i)%r(1:3)
3668 DO kk = 1, n_oxygens
3670 k = colvar%hydronium_dist_param%i_oxygens(kk)
3671 rpk(:) = my_particles(k)%r(1:3)
3672 rki =
pbc(rpk, rpi, cell)
3673 drki = sqrt(sum(rki**2))
3674 expfac_f_rki(ii) = expfac_f_rki(ii) + drki*expfac_f(kk)
3675 ddist_rki(1:3, kk, ii) = rki(1:3)/drki
3676 dexpfac_f_rki(kk, ii) = drki*dexpfac_f(kk)
3678 rion_num = rion_num + m(ii)*expfac_noh(ii)*expfac_f_rki(ii)
3682 rion_den = sum_expfac_noh*sum_expfac_f
3683 rion = rion_num/rion_den
3688 DO ii = 1, n_oxygens
3689 DO jj = 1, n_hydrogens
3690 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3691 + dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3692 *expfac_f_rki(ii)/rion_den
3693 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3694 - dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3695 *expfac_f_rki(ii)/rion_den
3696 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3697 + m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3698 *expfac_f_rki(ii)/rion_den
3699 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3700 - m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3701 *expfac_f_rki(ii)/rion_den
3703 DO kk = 1, n_oxygens
3705 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3706 - m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3707 *expfac_f(kk)/rion_den
3708 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3709 + m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3710 *expfac_f(kk)/rion_den
3711 DO jj = 1, n_hydrogens
3712 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3713 + m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3714 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3715 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3716 - m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3717 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3722 DO ii = 1, n_oxygens
3723 DO jj = 1, n_hydrogens
3724 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3725 - rion_num*sum_expfac_f*dexpfac_noh(ii) &
3726 *dnoh(1:3, jj, ii)/(rion_den**2)
3727 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3728 + rion_num*sum_expfac_f*dexpfac_noh(ii) &
3729 *dnoh(1:3, jj, ii)/(rion_den**2)
3730 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3731 - rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3732 *dnoh(1:3, jj, ii)/(rion_den**2)
3733 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3734 + rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3735 *dnoh(1:3, jj, ii)/(rion_den**2)
3739 DEALLOCATE (noh, m, f, expfac_noh, expfac_f)
3740 DEALLOCATE (dnoh, dm, df, dexpfac_noh, dexpfac_f)
3741 DEALLOCATE (ddist_rki, expfac_f_rki, dexpfac_f_rki)
3743 END SUBROUTINE hydronium_dist_colvar
3756 SUBROUTINE acid_hyd_dist_colvar(colvar, cell, subsys, particles)
3761 POINTER :: particles
3763 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3764 n_oxygens_acid, n_oxygens_water, &
3765 offseth, offseto, paoh, pcut, pwoh, &
3767 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac, expfac, nwoh
3768 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_rik
3769 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rik, dnaoh, dnwoh
3770 REAL(kind=
dp) :: dfcut, drik, drji, drjk, fbrace, fcut, fscalar, invden, invden_cut, lambda, &
3771 naoh, nc, num, num_cut, raoh, rik(3), rion, rion_den, rion_num, rji(3), rjk(3), rpi(3), &
3772 rpj(3), rpk(3), rrel, rwoh
3776 NULLIFY (my_particles, particles_i)
3778 n_oxygens_water = colvar%acid_hyd_dist_param%n_oxygens_water
3779 n_oxygens_acid = colvar%acid_hyd_dist_param%n_oxygens_acid
3780 n_hydrogens = colvar%acid_hyd_dist_param%n_hydrogens
3781 pwoh = colvar%acid_hyd_dist_param%pwoh
3782 qwoh = colvar%acid_hyd_dist_param%qwoh
3783 paoh = colvar%acid_hyd_dist_param%paoh
3784 qaoh = colvar%acid_hyd_dist_param%qaoh
3785 pcut = colvar%acid_hyd_dist_param%pcut
3786 qcut = colvar%acid_hyd_dist_param%qcut
3787 rwoh = colvar%acid_hyd_dist_param%rwoh
3788 raoh = colvar%acid_hyd_dist_param%raoh
3789 nc = colvar%acid_hyd_dist_param%nc
3790 lambda = colvar%acid_hyd_dist_param%lambda
3791 ALLOCATE (expfac(n_oxygens_water))
3792 ALLOCATE (nwoh(n_oxygens_water))
3793 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
3794 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
3795 ALLOCATE (dexpfac(n_oxygens_water))
3796 ALLOCATE (ddist_rik(3, n_oxygens_water, n_oxygens_acid))
3797 ALLOCATE (dexpfac_rik(n_oxygens_water, n_oxygens_acid))
3802 dnaoh(:, :, :) = 0._dp
3803 dnwoh(:, :, :) = 0._dp
3804 ddist_rik(:, :, :) = 0._dp
3806 dexpfac_rik(:, :) = 0._dp
3809 IF (
PRESENT(particles))
THEN
3810 my_particles => particles
3812 cpassert(
PRESENT(subsys))
3814 my_particles => particles_i%els
3818 DO ii = 1, n_oxygens_water
3819 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3820 rpi(:) = my_particles(i)%r(1:3)
3821 DO jj = 1, n_hydrogens
3822 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3823 rpj(:) = my_particles(j)%r(1:3)
3824 rji =
pbc(rpj, rpi, cell)
3825 drji = sqrt(sum(rji**2))
3827 num = 1.0_dp - rrel**pwoh
3828 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
3829 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3830 nwoh(ii) = nwoh(ii) + num*invden
3831 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
3832 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
3833 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3836 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
3837 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
3838 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3841 expfac(ii) = exp(lambda*nwoh(ii))
3842 dexpfac(ii) = lambda*expfac(ii)
3843 rion_den = rion_den + expfac(ii)
3847 DO kk = 1, n_oxygens_acid
3848 k = colvar%acid_hyd_dist_param%i_oxygens_acid(kk)
3849 rpk(:) = my_particles(k)%r(1:3)
3850 DO ii = 1, n_oxygens_water
3851 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3852 rpi(:) = my_particles(i)%r(1:3)
3853 rik =
pbc(rpi, rpk, cell)
3854 drik = sqrt(sum(rik**2))
3855 rion_num = rion_num + drik*expfac(ii)
3856 ddist_rik(1:3, ii, kk) = rik(1:3)/drik
3857 dexpfac_rik(ii, kk) = drik*dexpfac(ii)
3862 DO kk = 1, n_oxygens_acid
3863 k = colvar%acid_hyd_dist_param%i_oxygens_acid(kk)
3864 rpk(:) = my_particles(k)%r(1:3)
3865 DO jj = 1, n_hydrogens
3866 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3867 rpj(:) = my_particles(j)%r(1:3)
3868 rjk =
pbc(rpj, rpk, cell)
3869 drjk = sqrt(sum(rjk**2))
3871 num = 1.0_dp - rrel**paoh
3872 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
3873 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3874 naoh = naoh + num*invden
3875 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
3876 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
3877 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3880 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
3881 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
3882 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3886 num_cut = 1.0_dp - (naoh/nc)**pcut
3887 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
3888 fcut = num_cut*invden_cut
3892 fbrace = rion_num/rion_den/n_oxygens_acid
3897 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
3898 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
3899 offseto = n_oxygens_water
3900 offseth = n_oxygens_water + n_oxygens_acid
3901 DO kk = 1, n_oxygens_acid
3902 DO jj = 1, n_hydrogens
3903 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3904 + dfcut*dnaoh(1:3, jj, kk)*fbrace
3905 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3906 - dfcut*dnaoh(1:3, jj, kk)*fbrace
3912 DO kk = 1, n_oxygens_acid
3913 DO ii = 1, n_oxygens_water
3914 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3915 + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3917 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3918 - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3920 DO jj = 1, n_hydrogens
3921 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3922 + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3924 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3925 - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3931 DO ii = 1, n_oxygens_water
3932 DO jj = 1, n_hydrogens
3933 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3934 - fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3935 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3936 + fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3940 END SUBROUTINE acid_hyd_dist_colvar
3953 SUBROUTINE acid_hyd_shell_colvar(colvar, cell, subsys, particles)
3958 POINTER :: particles
3960 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, n_oxygens_acid, n_oxygens_water, offseth, &
3961 offseto, paoh, pcut, pm, poo, pwoh, qaoh, qcut, qm, qoo, qwoh, tt
3962 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dm, m, noo, nwoh, qloc
3963 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dnaoh, dnoo, dnwoh
3964 REAL(kind=
dp) :: dfcut, drji, drjk, drki, fcut, fscalar, invden, invden_cut, lambda, naoh, &
3965 nc, nh, num, num_cut, qsol, qtot, raoh, rji(3), rjk(3), rki(3), roo, rpi(3), rpj(3), &
3970 NULLIFY (my_particles, particles_i)
3972 n_oxygens_water = colvar%acid_hyd_shell_param%n_oxygens_water
3973 n_oxygens_acid = colvar%acid_hyd_shell_param%n_oxygens_acid
3974 n_hydrogens = colvar%acid_hyd_shell_param%n_hydrogens
3975 pwoh = colvar%acid_hyd_shell_param%pwoh
3976 qwoh = colvar%acid_hyd_shell_param%qwoh
3977 paoh = colvar%acid_hyd_shell_param%paoh
3978 qaoh = colvar%acid_hyd_shell_param%qaoh
3979 poo = colvar%acid_hyd_shell_param%poo
3980 qoo = colvar%acid_hyd_shell_param%qoo
3981 pm = colvar%acid_hyd_shell_param%pm
3982 qm = colvar%acid_hyd_shell_param%qm
3983 pcut = colvar%acid_hyd_shell_param%pcut
3984 qcut = colvar%acid_hyd_shell_param%qcut
3985 rwoh = colvar%acid_hyd_shell_param%rwoh
3986 raoh = colvar%acid_hyd_shell_param%raoh
3987 roo = colvar%acid_hyd_shell_param%roo
3988 nc = colvar%acid_hyd_shell_param%nc
3989 nh = colvar%acid_hyd_shell_param%nh
3990 lambda = colvar%acid_hyd_shell_param%lambda
3991 ALLOCATE (nwoh(n_oxygens_water))
3992 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
3993 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
3994 ALLOCATE (m(n_oxygens_water))
3995 ALLOCATE (dm(n_oxygens_water))
3996 ALLOCATE (noo(n_oxygens_water))
3997 ALLOCATE (dnoo(3, n_oxygens_water + n_oxygens_acid, n_oxygens_water))
3998 ALLOCATE (qloc(n_oxygens_water))
4002 dnaoh(:, :, :) = 0._dp
4003 dnwoh(:, :, :) = 0._dp
4004 dnoo(:, :, :) = 0._dp
4010 IF (
PRESENT(particles))
THEN
4011 my_particles => particles
4013 cpassert(
PRESENT(subsys))
4015 my_particles => particles_i%els
4019 DO ii = 1, n_oxygens_water
4020 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4021 rpi(:) = my_particles(i)%r(1:3)
4022 DO jj = 1, n_hydrogens
4023 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4024 rpj(:) = my_particles(j)%r(1:3)
4025 rji =
pbc(rpj, rpi, cell)
4026 drji = sqrt(sum(rji**2))
4028 num = 1.0_dp - rrel**pwoh
4029 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
4030 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4031 nwoh(ii) = nwoh(ii) + num*invden
4032 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
4033 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
4034 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4037 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
4038 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
4039 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4045 DO ii = 1, n_oxygens_water
4046 num = 1.0_dp - (nwoh(ii)/nh)**pm
4047 invden = 1.0_dp/(1.0_dp - (nwoh(ii)/nh)**qm)
4048 m(ii) = 1.0_dp - num*invden
4049 dm(ii) = (pm*(nwoh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
4050 (nwoh(ii)/nh)**(qm - 1))/nh
4054 DO ii = 1, n_oxygens_water
4055 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4056 rpi(:) = my_particles(i)%r(1:3)
4057 DO kk = 1, n_oxygens_water + n_oxygens_acid
4059 IF (kk <= n_oxygens_water)
THEN
4060 k = colvar%acid_hyd_shell_param%i_oxygens_water(kk)
4061 rpk(:) = my_particles(k)%r(1:3)
4063 tt = kk - n_oxygens_water
4064 k = colvar%acid_hyd_shell_param%i_oxygens_acid(tt)
4065 rpk(:) = my_particles(k)%r(1:3)
4067 rki =
pbc(rpk, rpi, cell)
4068 drki = sqrt(sum(rki**2))
4070 num = 1.0_dp - rrel**poo
4071 invden = 1.0_dp/(1.0_dp - rrel**qoo)
4072 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4073 noo(ii) = noo(ii) + num*invden
4074 fscalar = (-poo*(rrel**(poo - 1))*invden &
4075 + num*(invden**2)*qoo*(rrel**(qoo - 1)))/(drki*roo)
4076 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4079 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
4080 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drki)
4081 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4087 DO kk = 1, n_oxygens_acid
4088 k = colvar%acid_hyd_shell_param%i_oxygens_acid(kk)
4089 rpk(:) = my_particles(k)%r(1:3)
4090 DO jj = 1, n_hydrogens
4091 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4092 rpj(:) = my_particles(j)%r(1:3)
4093 rjk =
pbc(rpj, rpk, cell)
4094 drjk = sqrt(sum(rjk**2))
4096 num = 1.0_dp - rrel**paoh
4097 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
4098 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4099 naoh = naoh + num*invden
4100 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
4101 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
4102 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4105 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
4106 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
4107 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4111 num_cut = 1.0_dp - (naoh/nc)**pcut
4112 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
4113 fcut = num_cut*invden_cut
4116 DO ii = 1, n_oxygens_water
4117 qloc(ii) = exp(lambda*m(ii)*noo(ii))
4118 qtot = qtot + qloc(ii)
4120 qsol = log(qtot)/lambda
4121 colvar%ss = fcut*qsol
4124 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
4125 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
4126 offseto = n_oxygens_water
4127 offseth = n_oxygens_water + n_oxygens_acid
4128 DO kk = 1, n_oxygens_acid
4129 DO jj = 1, n_hydrogens
4130 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
4131 + dfcut*dnaoh(1:3, jj, kk)*qsol
4132 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4133 - dfcut*dnaoh(1:3, jj, kk)*qsol
4139 DO ii = 1, n_oxygens_water
4140 fscalar = fcut*qloc(ii)*dm(ii)*noo(ii)/qtot
4141 DO jj = 1, n_hydrogens
4142 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
4143 + fscalar*dnwoh(1:3, jj, ii)
4144 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4145 - fscalar*dnwoh(1:3, jj, ii)
4149 DO ii = 1, n_oxygens_water
4150 fscalar = fcut*qloc(ii)*m(ii)/qtot
4151 DO kk = 1, n_oxygens_water + n_oxygens_acid
4153 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + fscalar*dnoo(1:3, kk, ii)
4154 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) - fscalar*dnoo(1:3, kk, ii)
4158 END SUBROUTINE acid_hyd_shell_colvar
4170 SUBROUTINE coord_colvar(colvar, cell, subsys, particles)
4175 POINTER :: particles
4177 INTEGER :: i, ii, j, jj, k, kk, n_atoms_from, &
4178 n_atoms_to_a, n_atoms_to_b, p_a, p_b, &
4180 REAL(
dp) :: dfunc_ij, dfunc_jk, func_ij, func_jk, func_k, inv_n_atoms_from, invden_ij, &
4181 invden_jk, ncoord, num_ij, num_jk, r_0_a, r_0_b, rdist_ij, rdist_jk, rij, rjk
4182 REAL(
dp),
DIMENSION(3) :: ftmp_i, ftmp_j, ftmp_k, ss, xij, xjk, &
4190 NULLIFY (particles_i)
4192 IF (
PRESENT(particles))
THEN
4193 my_particles => particles
4195 cpassert(
PRESENT(subsys))
4197 my_particles => particles_i%els
4199 n_atoms_to_a = colvar%coord_param%n_atoms_to
4200 n_atoms_to_b = colvar%coord_param%n_atoms_to_b
4201 n_atoms_from = colvar%coord_param%n_atoms_from
4202 p_a = colvar%coord_param%nncrd
4203 q_a = colvar%coord_param%ndcrd
4204 r_0_a = colvar%coord_param%r_0
4205 p_b = colvar%coord_param%nncrd_b
4206 q_b = colvar%coord_param%ndcrd_b
4207 r_0_b = colvar%coord_param%r_0_b
4210 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
4211 DO ii = 1, n_atoms_from
4212 i = colvar%coord_param%i_at_from(ii)
4213 CALL get_coordinates(colvar, i, xpi, my_particles)
4214 DO jj = 1, n_atoms_to_a
4215 j = colvar%coord_param%i_at_to(jj)
4216 CALL get_coordinates(colvar, j, xpj, my_particles)
4219 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
4221 xij = matmul(cell%hmat, ss)
4222 rij = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
4223 IF (rij < 1.0e-8_dp) cycle
4224 rdist_ij = rij/r_0_a
4225 IF (abs(1.0_dp - rdist_ij) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4226 num_ij = (1.0_dp - rdist_ij**p_a)
4227 invden_ij = 1.0_dp/(1.0_dp - rdist_ij**q_a)
4228 func_ij = num_ij*invden_ij
4229 IF (rij < 1.0e-8_dp)
THEN
4233 dfunc_ij = (-p_a*rdist_ij**(p_a - 1)*invden_ij &
4234 + num_ij*(invden_ij)**2*q_a*rdist_ij**(q_a - 1))/(rij*r_0_a)
4238 func_ij = real(p_a, kind=
dp)/real(q_a, kind=
dp)
4239 dfunc_ij = real(p_a, kind=
dp)*real((-q_a + p_a), kind=
dp)/(real(2*q_a, kind=
dp)*r_0_a)
4241 IF (n_atoms_to_b /= 0)
THEN
4243 DO kk = 1, n_atoms_to_b
4244 k = colvar%coord_param%i_at_to_b(kk)
4246 CALL get_coordinates(colvar, k, xpk, my_particles)
4247 ss = matmul(cell%h_inv, xpj(:) - xpk(:))
4249 xjk = matmul(cell%hmat, ss)
4250 rjk = sqrt(xjk(1)**2 + xjk(2)**2 + xjk(3)**2)
4251 IF (rjk < 1.0e-8_dp) cycle
4252 rdist_jk = rjk/r_0_b
4253 IF (abs(1.0_dp - rdist_jk) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4254 num_jk = (1.0_dp - rdist_jk**p_b)
4255 invden_jk = 1.0_dp/(1.0_dp - rdist_jk**q_b)
4256 func_jk = num_jk*invden_jk
4257 IF (rjk < 1.0e-8_dp)
THEN
4261 dfunc_jk = (-p_b*rdist_jk**(p_b - 1)*invden_jk &
4262 + num_jk*(invden_jk)**2*q_b*rdist_jk**(q_b - 1))/(rjk*r_0_b)
4266 func_jk = real(p_b, kind=
dp)/real(q_b, kind=
dp)
4267 dfunc_jk = real(p_b, kind=
dp)*real((-q_b + p_b), kind=
dp)/(real(2*q_b, kind=
dp)*r_0_b)
4269 func_k = func_k + func_jk
4270 ftmp_k = -func_ij*dfunc_jk*xjk
4271 CALL put_derivative(colvar, n_atoms_from + n_atoms_to_a + kk, ftmp_k)
4273 ftmp_j = -dfunc_ij*xij*func_jk + func_ij*dfunc_jk*xjk
4274 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4279 ftmp_j = -dfunc_ij*xij
4280 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4282 ncoord = ncoord + func_ij*func_k
4283 ftmp_i = dfunc_ij*xij*func_k
4284 CALL put_derivative(colvar, ii, ftmp_i)
4287 colvar%ss = ncoord*inv_n_atoms_from
4288 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
4289 END SUBROUTINE coord_colvar
4298 SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
4304 POINTER :: particles
4306 INTEGER :: i, ii, j, jj, n_coord_from, n_coord_to, &
4308 REAL(
dp) :: den_n, den_q, fscalar, ftemp_i(3), inv_den_n, inv_den_q, lambda, num_n, num_q, &
4309 qfunc, r12, r_cut, rfact, rij(3), rpi(3), rpj(3)
4310 REAL(
dp),
DIMENSION(:),
POINTER :: dqfunc_dnl, expnl, nlcoord, sum_rij
4311 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: dnlcoord, dqfunc_dr
4318 NULLIFY (particles_i)
4320 IF (
PRESENT(particles))
THEN
4321 my_particles => particles
4323 cpassert(
PRESENT(subsys))
4325 my_particles => particles_i%els
4328 n_dist_from = colvar%mindist_param%n_dist_from
4329 n_coord_from = colvar%mindist_param%n_coord_from
4330 n_coord_to = colvar%mindist_param%n_coord_to
4331 p = colvar%mindist_param%p_exp
4332 q = colvar%mindist_param%q_exp
4333 r_cut = colvar%mindist_param%r_cut
4334 lambda = colvar%mindist_param%lambda
4336 NULLIFY (nlcoord, dnlcoord, dqfunc_dr, dqfunc_dnl, expnl, sum_rij)
4337 ALLOCATE (nlcoord(n_coord_from))
4338 ALLOCATE (dnlcoord(3, n_coord_from, n_coord_to))
4339 ALLOCATE (expnl(n_coord_from))
4340 ALLOCATE (sum_rij(n_coord_from))
4341 ALLOCATE (dqfunc_dr(3, n_dist_from, n_coord_from))
4342 ALLOCATE (dqfunc_dnl(n_coord_from))
4349 DO i = 1, n_coord_from
4350 ii = colvar%mindist_param%i_coord_from(i)
4351 rpi = my_particles(ii)%r(1:3)
4352 DO j = 1, n_coord_to
4353 jj = colvar%mindist_param%i_coord_to(j)
4354 rpj = my_particles(jj)%r(1:3)
4355 rij =
pbc(rpj, rpi, cell)
4356 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4358 num_n = 1.0_dp - rfact**p
4359 den_n = 1.0_dp - rfact**q
4360 inv_den_n = 1.0_dp/den_n
4361 IF (abs(inv_den_n) < 1.e-10_dp)
THEN
4362 inv_den_n = 1.e-10_dp
4366 fscalar = (-p*rfact**(p - 1) + num_n*q*rfact**(q - 1)*inv_den_n)*inv_den_n/(r_cut*r12)
4368 dnlcoord(1, i, j) = rij(1)*fscalar
4369 dnlcoord(2, i, j) = rij(2)*fscalar
4370 dnlcoord(3, i, j) = rij(3)*fscalar
4372 nlcoord(i) = nlcoord(i) + num_n*inv_den_n
4374 expnl(i) = exp(lambda*nlcoord(i))
4375 den_q = den_q + expnl(i)
4377 inv_den_q = 1.0_dp/den_q
4384 DO i = 1, n_dist_from
4385 ii = colvar%mindist_param%i_dist_from(i)
4386 rpi = my_particles(ii)%r(1:3)
4387 DO j = 1, n_coord_from
4388 jj = colvar%mindist_param%i_coord_from(j)
4389 rpj = my_particles(jj)%r(1:3)
4390 rij =
pbc(rpj, rpi, cell)
4391 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4393 num_q = num_q + r12*expnl(j)
4395 sum_rij(j) = sum_rij(j) + r12
4396 dqfunc_dr(1, i, j) = expnl(j)*rij(1)/r12
4397 dqfunc_dr(2, i, j) = expnl(j)*rij(2)/r12
4398 dqfunc_dr(3, i, j) = expnl(j)*rij(3)/r12
4405 qfunc = num_q*inv_den_q
4406 dqfunc_dr = dqfunc_dr*inv_den_q
4409 DO i = 1, n_coord_from
4410 dqfunc_dnl(i) = lambda*expnl(i)*inv_den_q*(sum_rij(i) - num_q*inv_den_q)
4414 DO i = 1, n_dist_from
4415 DO j = 1, n_coord_from
4416 ftemp_i(1) = dqfunc_dr(1, i, j)
4417 ftemp_i(2) = dqfunc_dr(2, i, j)
4418 ftemp_i(3) = dqfunc_dr(3, i, j)
4420 CALL put_derivative(colvar, i, ftemp_i)
4421 CALL put_derivative(colvar, j + n_dist_from, -ftemp_i)
4425 DO i = 1, n_coord_from
4426 DO j = 1, n_coord_to
4427 ftemp_i(1) = dqfunc_dnl(i)*dnlcoord(1, i, j)
4428 ftemp_i(2) = dqfunc_dnl(i)*dnlcoord(2, i, j)
4429 ftemp_i(3) = dqfunc_dnl(i)*dnlcoord(3, i, j)
4431 CALL put_derivative(colvar, i + n_dist_from, ftemp_i)
4432 CALL put_derivative(colvar, j + n_dist_from + n_coord_from, -ftemp_i)
4437 DEALLOCATE (nlcoord)
4438 DEALLOCATE (dnlcoord)
4440 DEALLOCATE (dqfunc_dr)
4441 DEALLOCATE (sum_rij)
4442 DEALLOCATE (dqfunc_dnl)
4444 END SUBROUTINE mindist_colvar
4454 SUBROUTINE combine_colvar(colvar, cell, subsys, particles)
4459 POINTER :: particles
4461 CHARACTER(LEN=default_string_length) :: def_error, this_error
4462 CHARACTER(LEN=default_string_length), &
4463 ALLOCATABLE,
DIMENSION(:) :: my_par
4464 INTEGER :: i, ii, j, ncolv, ndim
4466 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dss_vals, my_val, ss_vals
4467 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi
4472 IF (
PRESENT(particles))
THEN
4473 my_particles => particles
4475 cpassert(
PRESENT(subsys))
4477 my_particles => particles_i%els
4480 ncolv =
SIZE(colvar%combine_cvs_param%colvar_p)
4481 ALLOCATE (ss_vals(ncolv))
4482 ALLOCATE (dss_vals(ncolv))
4486 CALL colvar_recursive_eval(colvar%combine_cvs_param%colvar_p(i)%colvar, cell, my_particles)
4487 ss_vals(i) = colvar%combine_cvs_param%colvar_p(i)%colvar%ss
4492 ndim =
SIZE(colvar%combine_cvs_param%c_parameters) + &
4493 SIZE(colvar%combine_cvs_param%variables)
4494 ALLOCATE (my_par(ndim))
4495 my_par(1:
SIZE(colvar%combine_cvs_param%variables)) = colvar%combine_cvs_param%variables
4496 my_par(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%c_parameters
4497 ALLOCATE (my_val(ndim))
4498 my_val(1:
SIZE(colvar%combine_cvs_param%variables)) = ss_vals
4499 my_val(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%v_parameters
4500 CALL parsef(1, trim(colvar%combine_cvs_param%function), my_par)
4501 colvar%ss =
evalf(1, my_val)
4503 dss_vals(i) =
evalfd(1, i, my_val, colvar%combine_cvs_param%dx, err)
4504 IF ((abs(err) > colvar%combine_cvs_param%lerr))
THEN
4505 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
4506 WRITE (def_error,
"(A,G12.6,A)")
"(", colvar%combine_cvs_param%lerr,
")"
4509 CALL cp_warn(__location__, &
4510 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
4511 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
4512 trim(def_error)//
' . ')
4520 ALLOCATE (fi(3, colvar%n_atom_s))
4523 DO j = 1, colvar%combine_cvs_param%colvar_p(i)%colvar%n_atom_s
4525 fi(:, ii) = colvar%combine_cvs_param%colvar_p(i)%colvar%dsdr(:, j)*dss_vals(i)
4529 DO i = 1, colvar%n_atom_s
4530 CALL put_derivative(colvar, i, fi(:, i))
4534 DEALLOCATE (ss_vals)
4535 DEALLOCATE (dss_vals)
4536 END SUBROUTINE combine_colvar
4552 SUBROUTINE reaction_path_colvar(colvar, cell, subsys, particles)
4557 POINTER :: particles
4563 IF (
PRESENT(particles))
THEN
4564 my_particles => particles
4566 cpassert(
PRESENT(subsys))
4568 my_particles => particles_i%els
4571 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4572 CALL rpath_dist_rmsd(colvar, my_particles)
4573 ELSE IF (colvar%reaction_path_param%rmsd)
THEN
4574 CALL rpath_rmsd(colvar, my_particles)
4576 CALL rpath_colvar(colvar, cell, my_particles)
4579 END SUBROUTINE reaction_path_colvar
4590 SUBROUTINE rpath_colvar(colvar, cell, particles)
4595 INTEGER :: i, iend, ii, istart, j, k, ncolv, nconf
4596 REAL(
dp) :: lambda, step_size
4597 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: s1, ss_vals
4598 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, f_vals, fi, s1v
4599 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
4601 istart = colvar%reaction_path_param%function_bounds(1)
4602 iend = colvar%reaction_path_param%function_bounds(2)
4604 nconf = colvar%reaction_path_param%nr_frames
4605 step_size = colvar%reaction_path_param%step_size
4606 ncolv = colvar%reaction_path_param%n_components
4607 lambda = colvar%reaction_path_param%lambda
4608 ALLOCATE (f_vals(ncolv, istart:iend))
4609 f_vals(:, :) = colvar%reaction_path_param%f_vals
4610 ALLOCATE (ss_vals(ncolv))
4613 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4614 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4617 ALLOCATE (s1v(2, istart:iend))
4618 ALLOCATE (ds1v(ncolv, 2, istart:iend))
4621 ALLOCATE (ds1(ncolv, 2))
4624 s1v(1, k) = real(k, kind=
dp)*step_size*exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4625 s1v(2, k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4627 ds1v(j, 1, k) = f_vals(j, k)*s1v(1, k)
4628 ds1v(j, 2, k) = f_vals(j, k)*s1v(2, k)
4638 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4640 ALLOCATE (fi(3, colvar%n_atom_s))
4644 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
4646 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)*lambda* &
4647 (ds1(i, 1)/s1(2)/real(nconf - 1,
dp) - colvar%ss*ds1(i, 2)/s1(2))*2.0_dp
4651 DO i = 1, colvar%n_atom_s
4652 CALL put_derivative(colvar, i, fi(:, i))
4657 DEALLOCATE (ss_vals)
4663 END SUBROUTINE rpath_colvar
4674 SUBROUTINE rpath_dist_rmsd(colvar, particles)
4678 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4679 INTEGER,
DIMENSION(:),
POINTER :: iatom
4680 REAL(
dp) :: lambda, my_rmsd, s1(2), sum_exp
4681 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, vec_dif
4682 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dvec_dif, fi, riat, s1v
4683 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4684 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4685 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4687 nconf = colvar%reaction_path_param%nr_frames
4688 rmsd_atom = colvar%reaction_path_param%n_components
4689 lambda = colvar%reaction_path_param%lambda
4690 path_conf => colvar%reaction_path_param%r_ref
4691 iatom => colvar%reaction_path_param%i_rmsd
4693 natom =
SIZE(particles)
4695 ALLOCATE (r0(3*natom))
4696 ALLOCATE (r(3*natom))
4697 ALLOCATE (riat(3, rmsd_atom))
4698 ALLOCATE (vec_dif(rmsd_atom))
4699 ALLOCATE (dvec_dif(3, rmsd_atom))
4700 ALLOCATE (s1v(2, nconf))
4701 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4702 ALLOCATE (ds1(3, rmsd_atom, 2))
4705 r0(ii + 1) = particles(i)%r(1)
4706 r0(ii + 2) = particles(i)%r(2)
4707 r0(ii + 3) = particles(i)%r(3)
4710 DO iat = 1, rmsd_atom
4712 riat(:, iat) = particles(ii)%r
4718 r(ii + 1) = path_conf(ii + 1, ik)
4719 r(ii + 2) = path_conf(ii + 2, ik)
4720 r(ii + 3) = path_conf(ii + 3, ik)
4723 CALL rmsd3(particles, r, r0, output_unit=-1, my_val=my_rmsd, rotate=.true.)
4726 DO iat = 1, rmsd_atom
4729 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 &
4730 + (riat(3, iat) - r(ii + 3))**2
4731 sum_exp = sum_exp + vec_dif(iat)
4734 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*sum_exp)
4735 s1v(2, ik) = exp(-lambda*sum_exp)
4736 DO iat = 1, rmsd_atom
4739 ds1v(1, iat, 1, ik) = r(ii + 1)*s1v(1, ik)
4740 ds1v(1, iat, 2, ik) = r(ii + 1)*s1v(2, ik)
4741 ds1v(2, iat, 1, ik) = r(ii + 2)*s1v(1, ik)
4742 ds1v(2, iat, 2, ik) = r(ii + 2)*s1v(2, ik)
4743 ds1v(3, iat, 1, ik) = r(ii + 3)*s1v(1, ik)
4744 ds1v(3, iat, 2, ik) = r(ii + 3)*s1v(2, ik)
4751 DO iat = 1, rmsd_atom
4758 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4760 ALLOCATE (fi(3, rmsd_atom))
4762 DO iat = 1, rmsd_atom
4763 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))
4764 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))
4765 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))
4766 CALL put_derivative(colvar, iat, fi(:, iat))
4773 DEALLOCATE (vec_dif)
4774 DEALLOCATE (dvec_dif)
4779 END SUBROUTINE rpath_dist_rmsd
4786 SUBROUTINE rpath_rmsd(colvar, particles)
4790 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4791 INTEGER,
DIMENSION(:),
POINTER :: iatom
4792 REAL(
dp) :: lambda, my_rmsd, s1(2)
4793 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0
4794 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi, riat, s1v
4795 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4796 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4797 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4798 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
4799 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
4801 nconf = colvar%reaction_path_param%nr_frames
4802 rmsd_atom = colvar%reaction_path_param%n_components
4803 lambda = colvar%reaction_path_param%lambda
4804 path_conf => colvar%reaction_path_param%r_ref
4805 iatom => colvar%reaction_path_param%i_rmsd
4807 natom =
SIZE(particles)
4809 ALLOCATE (r0(3*natom))
4810 ALLOCATE (r(3*natom))
4811 ALLOCATE (riat(3, rmsd_atom))
4812 ALLOCATE (s1v(2, nconf))
4813 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4814 ALLOCATE (ds1(3, rmsd_atom, 2))
4815 ALLOCATE (drmsd(3, natom))
4817 ALLOCATE (weight(natom))
4821 r0(ii + 1) = particles(i)%r(1)
4822 r0(ii + 2) = particles(i)%r(2)
4823 r0(ii + 3) = particles(i)%r(3)
4826 DO iat = 1, rmsd_atom
4828 riat(:, iat) = particles(ii)%r
4833 DO iat = 1, rmsd_atom
4841 r(ii + 1) = path_conf(ii + 1, ik)
4842 r(ii + 2) = path_conf(ii + 2, ik)
4843 r(ii + 3) = path_conf(ii + 3, ik)
4846 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
4847 rotate=.false., drmsd3=drmsd)
4849 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*my_rmsd)
4850 s1v(2, ik) = exp(-lambda*my_rmsd)
4851 DO iat = 1, rmsd_atom
4853 ds1v(1, iat, 1, ik) = drmsd(1, i)*s1v(1, ik)
4854 ds1v(1, iat, 2, ik) = drmsd(1, i)*s1v(2, ik)
4855 ds1v(2, iat, 1, ik) = drmsd(2, i)*s1v(1, ik)
4856 ds1v(2, iat, 2, ik) = drmsd(2, i)*s1v(2, ik)
4857 ds1v(3, iat, 1, ik) = drmsd(3, i)*s1v(1, ik)
4858 ds1v(3, iat, 2, ik) = drmsd(3, i)*s1v(2, ik)
4865 DO iat = 1, rmsd_atom
4872 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4874 ALLOCATE (fi(3, rmsd_atom))
4876 DO iat = 1, rmsd_atom
4877 fi(1, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(1, iat, 1) - ds1(1, iat, 2)*s1(1)/s1(2))
4878 fi(2, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(2, iat, 1) - ds1(2, iat, 2)*s1(1)/s1(2))
4879 fi(3, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(3, iat, 1) - ds1(3, iat, 2)*s1(1)/s1(2))
4880 CALL put_derivative(colvar, iat, fi(:, iat))
4893 END SUBROUTINE rpath_rmsd
4905 SUBROUTINE distance_from_path_colvar(colvar, cell, subsys, particles)
4910 POINTER :: particles
4916 IF (
PRESENT(particles))
THEN
4917 my_particles => particles
4919 cpassert(
PRESENT(subsys))
4921 my_particles => particles_i%els
4924 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4925 CALL dpath_dist_rmsd(colvar, my_particles)
4926 ELSE IF (colvar%reaction_path_param%rmsd)
THEN
4927 CALL dpath_rmsd(colvar, my_particles)
4929 CALL dpath_colvar(colvar, cell, my_particles)
4932 END SUBROUTINE distance_from_path_colvar
4944 SUBROUTINE dpath_colvar(colvar, cell, particles)
4949 INTEGER :: i, iend, ii, istart, j, k, ncolv
4950 REAL(
dp) :: lambda, s1
4951 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: ds1, s1v, ss_vals
4952 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1v, f_vals, fi
4954 istart = colvar%reaction_path_param%function_bounds(1)
4955 iend = colvar%reaction_path_param%function_bounds(2)
4957 ncolv = colvar%reaction_path_param%n_components
4958 lambda = colvar%reaction_path_param%lambda
4959 ALLOCATE (f_vals(ncolv, istart:iend))
4960 f_vals(:, :) = colvar%reaction_path_param%f_vals
4961 ALLOCATE (ss_vals(ncolv))
4964 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4965 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4968 ALLOCATE (s1v(istart:iend))
4969 ALLOCATE (ds1v(ncolv, istart:iend))
4970 ALLOCATE (ds1(ncolv))
4973 s1v(k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4975 ds1v(j, k) = f_vals(j, k)*s1v(k)
4983 colvar%ss = -1.0_dp/lambda*log(s1)
4985 ALLOCATE (fi(3, colvar%n_atom_s))
4989 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
4991 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)* &
4992 2.0_dp*(ss_vals(i) - ds1(i)/s1)
4996 DO i = 1, colvar%n_atom_s
4997 CALL put_derivative(colvar, i, fi(:, i))
5002 DEALLOCATE (ss_vals)
5007 END SUBROUTINE dpath_colvar
5018 SUBROUTINE dpath_dist_rmsd(colvar, particles)
5023 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5024 INTEGER,
DIMENSION(:),
POINTER :: iatom
5025 REAL(
dp) :: lambda, s1, sum_exp
5026 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v, vec_dif
5027 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, dvec_dif, fi, riat
5028 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5029 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5031 nconf = colvar%reaction_path_param%nr_frames
5032 rmsd_atom = colvar%reaction_path_param%n_components
5033 lambda = colvar%reaction_path_param%lambda
5034 path_conf => colvar%reaction_path_param%r_ref
5035 iatom => colvar%reaction_path_param%i_rmsd
5037 natom =
SIZE(particles)
5039 ALLOCATE (r0(3*natom))
5040 ALLOCATE (r(3*natom))
5041 ALLOCATE (riat(3, rmsd_atom))
5042 ALLOCATE (vec_dif(rmsd_atom))
5043 ALLOCATE (dvec_dif(3, rmsd_atom))
5044 ALLOCATE (s1v(nconf))
5045 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5046 ALLOCATE (ds1(3, rmsd_atom))
5049 r0(ii + 1) = particles(i)%r(1)
5050 r0(ii + 2) = particles(i)%r(2)
5051 r0(ii + 3) = particles(i)%r(3)
5054 DO iat = 1, rmsd_atom
5056 riat(:, iat) = particles(ii)%r
5062 r(ii + 1) = path_conf(ii + 1, ik)
5063 r(ii + 2) = path_conf(ii + 2, ik)
5064 r(ii + 3) = path_conf(ii + 3, ik)
5067 CALL rmsd3(particles, r, r0, output_unit=-1, rotate=.true.)
5070 DO iat = 1, rmsd_atom
5073 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 + (riat(3, iat) - r(ii + 3))**2
5074 sum_exp = sum_exp + vec_dif(iat)
5075 dvec_dif(1, iat) = r(ii + 1)
5076 dvec_dif(2, iat) = r(ii + 2)
5077 dvec_dif(3, iat) = r(ii + 3)
5079 s1v(ik) = exp(-lambda*sum_exp)
5080 DO iat = 1, rmsd_atom
5081 ds1v(1, iat, ik) = dvec_dif(1, iat)*s1v(ik)
5082 ds1v(2, iat, ik) = dvec_dif(2, iat)*s1v(ik)
5083 ds1v(3, iat, ik) = dvec_dif(3, iat)*s1v(ik)
5088 DO iat = 1, rmsd_atom
5093 colvar%ss = -1.0_dp/lambda*log(s1)
5095 ALLOCATE (fi(3, rmsd_atom))
5097 DO iat = 1, rmsd_atom
5098 fi(:, iat) = 2.0_dp*(riat(:, iat) - ds1(:, iat)/s1)
5099 CALL put_derivative(colvar, iat, fi(:, iat))
5106 DEALLOCATE (vec_dif)
5107 DEALLOCATE (dvec_dif)
5111 END SUBROUTINE dpath_dist_rmsd
5118 SUBROUTINE dpath_rmsd(colvar, particles)
5123 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5124 INTEGER,
DIMENSION(:),
POINTER :: iatom
5125 REAL(
dp) :: lambda, my_rmsd, s1
5126 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v
5127 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, fi, riat
5128 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5129 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5130 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
5131 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
5133 nconf = colvar%reaction_path_param%nr_frames
5134 rmsd_atom = colvar%reaction_path_param%n_components
5135 lambda = colvar%reaction_path_param%lambda
5136 path_conf => colvar%reaction_path_param%r_ref
5137 iatom => colvar%reaction_path_param%i_rmsd
5139 natom =
SIZE(particles)
5141 ALLOCATE (r0(3*natom))
5142 ALLOCATE (r(3*natom))
5143 ALLOCATE (riat(3, rmsd_atom))
5144 ALLOCATE (s1v(nconf))
5145 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5146 ALLOCATE (ds1(3, rmsd_atom))
5147 ALLOCATE (drmsd(3, natom))
5149 ALLOCATE (weight(natom))
5153 r0(ii + 1) = particles(i)%r(1)
5154 r0(ii + 2) = particles(i)%r(2)
5155 r0(ii + 3) = particles(i)%r(3)
5158 DO iat = 1, rmsd_atom
5160 riat(:, iat) = particles(ii)%r
5165 DO iat = 1, rmsd_atom
5173 r(ii + 1) = path_conf(ii + 1, ik)
5174 r(ii + 2) = path_conf(ii + 2, ik)
5175 r(ii + 3) = path_conf(ii + 3, ik)
5178 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
5179 rotate=.false., drmsd3=drmsd)
5181 s1v(ik) = exp(-lambda*my_rmsd)
5182 DO iat = 1, rmsd_atom
5184 ds1v(1, iat, ik) = drmsd(1, i)*s1v(ik)
5185 ds1v(2, iat, ik) = drmsd(2, i)*s1v(ik)
5186 ds1v(3, iat, ik) = drmsd(3, i)*s1v(ik)
5191 DO iat = 1, rmsd_atom
5196 colvar%ss = -1.0_dp/lambda*log(s1)
5198 ALLOCATE (fi(3, rmsd_atom))
5200 DO iat = 1, rmsd_atom
5201 fi(:, iat) = ds1(:, iat)/s1
5202 CALL put_derivative(colvar, iat, fi(:, iat))
5215 END SUBROUTINE dpath_rmsd
5226 SUBROUTINE population_colvar(colvar, cell, subsys, particles)
5231 POINTER :: particles
5233 INTEGER :: i, ii, jj, n_atoms_from, n_atoms_to, &
5235 REAL(
dp) :: dfunc, dfunc_coord, ftmp(3), func, func_coord, inv_n_atoms_from, invden, n_0, &
5236 ncoord, norm, num, population, r12, r_0, rdist, sigma, ss(3), xij(3)
5237 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftmp_coord
5238 REAL(
dp),
DIMENSION(3) :: xpi, xpj
5245 NULLIFY (particles_i)
5247 IF (
PRESENT(particles))
THEN
5248 my_particles => particles
5250 cpassert(
PRESENT(subsys))
5252 my_particles => particles_i%els
5254 n_atoms_to = colvar%population_param%n_atoms_to
5255 n_atoms_from = colvar%population_param%n_atoms_from
5256 nncrd = colvar%population_param%nncrd
5257 ndcrd = colvar%population_param%ndcrd
5258 r_0 = colvar%population_param%r_0
5259 n_0 = colvar%population_param%n0
5260 sigma = colvar%population_param%sigma
5262 ALLOCATE (ftmp_coord(3, n_atoms_to))
5268 colvar%dsdr = 0.0_dp
5269 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
5271 norm = sqrt(
pi*2.0_dp)*sigma
5274 DO ii = 1, n_atoms_from
5275 i = colvar%population_param%i_at_from(ii)
5276 CALL get_coordinates(colvar, i, xpi, my_particles)
5277 DO jj = 1, n_atoms_to
5278 i = colvar%population_param%i_at_to(jj)
5279 CALL get_coordinates(colvar, i, xpj, my_particles)
5280 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
5282 xij = matmul(cell%hmat, ss)
5283 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
5284 IF (r12 < 1.0e-8_dp) cycle
5286 num = (1.0_dp - rdist**nncrd)
5287 invden = 1.0_dp/(1.0_dp - rdist**ndcrd)
5288 func_coord = num*invden
5289 dfunc_coord = (-nncrd*rdist**(nncrd - 1)*invden &
5290 + num*(invden)**2*ndcrd*rdist**(ndcrd - 1))/(r12*r_0)
5292 ncoord = ncoord + func_coord
5293 ftmp_coord(1, jj) = dfunc_coord*xij(1)
5294 ftmp_coord(2, jj) = dfunc_coord*xij(2)
5295 ftmp_coord(3, jj) = dfunc_coord*xij(3)
5298 func = exp(-(ncoord - n_0)**2/(2.0_dp*sigma*sigma))
5299 dfunc = -func*(ncoord - n_0)/(sigma*sigma)
5301 population = population + norm*func
5302 DO jj = 1, n_atoms_to
5303 ftmp(1) = ftmp_coord(1, jj)*dfunc
5304 ftmp(2) = ftmp_coord(2, jj)*dfunc
5305 ftmp(3) = ftmp_coord(3, jj)*dfunc
5306 CALL put_derivative(colvar, ii, ftmp)
5307 ftmp(1) = -ftmp_coord(1, jj)*dfunc
5308 ftmp(2) = -ftmp_coord(2, jj)*dfunc
5309 ftmp(3) = -ftmp_coord(3, jj)*dfunc
5310 CALL put_derivative(colvar, n_atoms_from + jj, ftmp)
5314 colvar%ss = population
5315 END SUBROUTINE population_colvar
5327 SUBROUTINE gyration_radius_colvar(colvar, cell, subsys, particles)
5333 POINTER :: particles
5335 INTEGER :: i, ii, n_atoms
5336 REAL(
dp) :: dri2, func, gyration, inv_n, mass_tot, mi
5337 REAL(
dp),
DIMENSION(3) :: dfunc, dxi, ftmp, ss, xpcom, xpi
5341 NULLIFY (particles_i, my_particles)
5343 IF (
PRESENT(particles))
THEN
5344 my_particles => particles
5346 cpassert(
PRESENT(subsys))
5348 my_particles => particles_i%els
5350 n_atoms = colvar%gyration_param%n_atoms
5351 inv_n = 1.0_dp/n_atoms
5357 i = colvar%gyration_param%i_at(ii)
5358 CALL get_coordinates(colvar, i, xpi, my_particles)
5359 CALL get_mass(colvar, i, mi, my_particles)
5360 xpcom(:) = xpcom(:) + xpi(:)*mi
5361 mass_tot = mass_tot + mi
5363 xpcom(:) = xpcom(:)/mass_tot
5369 i = colvar%gyration_param%i_at(ii)
5370 CALL get_coordinates(colvar, i, xpi, my_particles)
5371 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5373 dxi = matmul(cell%hmat, ss)
5374 dri2 = (dxi(1)**2 + dxi(2)**2 + dxi(3)**2)
5376 dfunc(:) = dfunc(:) + dxi(:)
5378 gyration = sqrt(inv_n*func)
5381 i = colvar%gyration_param%i_at(ii)
5382 CALL get_coordinates(colvar, i, xpi, my_particles)
5383 CALL get_mass(colvar, i, mi, my_particles)
5384 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5386 dxi = matmul(cell%hmat, ss)
5387 ftmp(1) = dxi(1) - dfunc(1)*mi/mass_tot
5388 ftmp(2) = dxi(2) - dfunc(2)*mi/mass_tot
5389 ftmp(3) = dxi(3) - dfunc(3)*mi/mass_tot
5390 ftmp(:) = ftmp(:)*inv_n/gyration
5391 CALL put_derivative(colvar, ii, ftmp)
5393 colvar%ss = gyration
5395 END SUBROUTINE gyration_radius_colvar
5406 SUBROUTINE rmsd_colvar(colvar, subsys, particles)
5410 POINTER :: particles
5412 CALL rmsd_colvar_low(colvar, subsys, particles)
5413 END SUBROUTINE rmsd_colvar
5428 SUBROUTINE rmsd_colvar_low(colvar, subsys, particles)
5433 POINTER :: particles
5435 INTEGER :: i, ii, natom, nframes
5436 REAL(kind=
dp) :: cv_val, f1, ftmp(3)
5437 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: der, r,
rmsd
5438 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r0
5439 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: drmsd
5440 REAL(kind=
dp),
DIMENSION(:),
POINTER :: weights
5444 NULLIFY (my_particles, particles_i, weights)
5446 IF (
PRESENT(particles))
THEN
5447 my_particles => particles
5449 cpassert(
PRESENT(subsys))
5451 my_particles => particles_i%els
5454 natom =
SIZE(my_particles)
5455 nframes = colvar%rmsd_param%nr_frames
5456 ALLOCATE (drmsd(3, natom, nframes))
5459 ALLOCATE (r0(3*natom, nframes))
5460 ALLOCATE (
rmsd(nframes))
5461 ALLOCATE (der(nframes))
5462 ALLOCATE (r(3*natom))
5464 weights => colvar%rmsd_param%weights
5467 r(ii + 1) = my_particles(i)%r(1)
5468 r(ii + 2) = my_particles(i)%r(2)
5469 r(ii + 3) = my_particles(i)%r(3)
5471 r0(:, :) = colvar%rmsd_param%r_ref
5474 CALL rmsd3(my_particles, r, r0(:, 1), output_unit=-1, weights=weights, my_val=
rmsd(1), rotate=.false., drmsd3=drmsd(:, :, 1))
5476 IF (nframes == 2)
THEN
5477 CALL rmsd3(my_particles, r, r0(:, 2), output_unit=-1, weights=weights, &
5478 my_val=
rmsd(2), rotate=.false., drmsd3=drmsd(:, :, 2))
5484 der(1) = f1 - cv_val*f1
5486 der(2) = -f1 - cv_val*f1
5488 DO i = 1, colvar%rmsd_param%n_atoms
5489 ii = colvar%rmsd_param%i_rmsd(i)
5490 IF (weights(ii) > 0.0_dp)
THEN
5491 ftmp(1) = der(1)*drmsd(1, ii, 1) + der(2)*drmsd(1, ii, 2)
5492 ftmp(2) = der(1)*drmsd(2, ii, 1) + der(2)*drmsd(2, ii, 2)
5493 ftmp(3) = der(1)*drmsd(3, ii, 1) + der(2)*drmsd(3, ii, 2)
5494 CALL put_derivative(colvar, i, ftmp)
5497 ELSE IF (nframes == 1)
THEN
5500 cv_val = sqrt(
rmsd(1))
5502 IF (cv_val /= 0.0_dp) f1 = 0.5_dp/cv_val
5503 DO i = 1, colvar%rmsd_param%n_atoms
5504 ii = colvar%rmsd_param%i_rmsd(i)
5505 IF (weights(ii) > 0.0_dp)
THEN
5506 ftmp(1) = f1*drmsd(1, ii, 1)
5507 ftmp(2) = f1*drmsd(2, ii, 1)
5508 ftmp(3) = f1*drmsd(3, ii, 1)
5509 CALL put_derivative(colvar, i, ftmp)
5513 cpabort(
"RMSD implemented only for 1 and 2 reference frames!")
5523 END SUBROUTINE rmsd_colvar_low
5535 SUBROUTINE ring_puckering_colvar(colvar, cell, subsys, particles)
5540 POINTER :: particles
5542 INTEGER :: i, ii, j, jj, m, nring
5543 REAL(kind=
dp) :: a, at, b, da, db, ds, kr, rpxpp, svar
5544 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cosj, sinj, z
5545 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r
5546 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: nforce, zforce
5547 REAL(kind=
dp),
DIMENSION(3) :: ftmp, nv, r0, rp, rpp, uv
5548 REAL(kind=
dp),
DIMENSION(3, 3) :: dnvp, dnvpp
5553 IF (
PRESENT(particles))
THEN
5554 my_particles => particles
5556 cpassert(
PRESENT(subsys))
5558 my_particles => particles_i%els
5561 nring = colvar%ring_puckering_param%nring
5562 ALLOCATE (r(3, nring), z(nring), cosj(nring), sinj(nring))
5563 ALLOCATE (nforce(3, 3, nring), zforce(nring, nring, 3))
5565 i = colvar%ring_puckering_param%atoms(ii)
5566 CALL get_coordinates(colvar, i, r(:, ii), my_particles)
5571 r(:, ii) =
pbc(r(:, ii), r0, cell)
5576 r0(:) = r0(:) + r(:, ii)
5578 kr = 1._dp/real(nring, kind=
dp)
5581 r(:, ii) = r(:, ii) - r0(:)
5587 cosj(ii) = cos(
twopi*(ii - 1)*kr)
5588 sinj(ii) = sin(
twopi*(ii - 1)*kr)
5589 rp(:) = rp(:) + r(:, ii)*sinj(ii)
5590 rpp(:) = rpp(:) + r(:, ii)*cosj(ii)
5593 nv = nv/sqrt(sum(nv**2))
5597 rpxpp = sqrt(sum(uv**2))
5602 dnvp(:, i) = uv - nv*sum(uv*nv)
5606 dnvpp(:, i) = uv - nv*sum(uv*nv)
5609 nforce(:, :, ii) = dnvp(:, :)*sinj(ii) + dnvpp(:, :)*cosj(ii)
5614 z(ii) = sum(r(:, ii)*nv(:))
5620 zforce(ii, jj, :) = nv
5622 zforce(ii, jj, :) = 0._dp
5626 zforce(ii, jj, i) = zforce(ii, jj, i) + r(j, ii)*nforce(j, i, jj)
5632 IF (colvar%ring_puckering_param%iq == 0)
THEN
5634 svar = sqrt(sum(z**2))
5638 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*z(jj)
5641 CALL put_derivative(colvar, ii, ftmp)
5644 m = abs(colvar%ring_puckering_param%iq)
5646 IF (mod(nring, 2) == 0 .AND. colvar%ring_puckering_param%iq == nring/2)
THEN
5650 IF (mod(ii, 2) == 0)
THEN
5656 svar = svar*sqrt(kr)
5660 IF (mod(jj, 2) == 0)
THEN
5661 ftmp(:) = ftmp(:) - zforce(jj, ii, :)*sqrt(kr)
5663 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*sqrt(kr)
5666 CALL put_derivative(colvar, ii, -ftmp)
5669 cpassert(m <= (nring - 1)/2)
5673 a = a + z(ii)*cos(
twopi*m*(ii - 1)*kr)
5674 b = b - z(ii)*sin(
twopi*m*(ii - 1)*kr)
5676 a = a*sqrt(2._dp*kr)
5677 b = b*sqrt(2._dp*kr)
5678 IF (colvar%ring_puckering_param%iq > 0)
THEN
5680 svar = sqrt(a*a + b*b)
5686 IF (at >
pi/2._dp)
THEN
5687 svar = 2.5_dp*
pi - at
5689 svar = 0.5_dp*
pi - at
5697 ds = da*cos(
twopi*m*(ii - 1)*kr)
5698 ds = ds - db*sin(
twopi*m*(ii - 1)*kr)
5699 ftmp(:) = ftmp(:) + ds*sqrt(2._dp*kr)*zforce(ii, jj, :)
5701 CALL put_derivative(colvar, jj, ftmp)
5708 DEALLOCATE (r, z, cosj, sinj, nforce, zforce)
5710 END SUBROUTINE ring_puckering_colvar
5732 RECURSIVE FUNCTION rec_eval_grid(iw1, ncol, f_vals, v_count, &
5733 gp, grid_sp, step_size, istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
RESULT(k)
5734 INTEGER :: iw1, ncol
5735 REAL(
dp),
DIMENSION(:, :),
POINTER :: f_vals
5737 REAL(
dp),
DIMENSION(:),
POINTER :: gp, grid_sp
5738 REAL(
dp) :: step_size
5739 INTEGER :: istart, iend
5740 REAL(
dp),
DIMENSION(:, :),
POINTER :: s1v
5741 REAL(
dp),
DIMENSION(:),
POINTER :: s1
5742 INTEGER,
DIMENSION(:, :),
POINTER :: p_bounds
5744 INTEGER :: ifunc, nconf, k
5746 INTEGER :: count1, i
5749 IF (v_count < ncol)
THEN
5750 count1 = v_count + 1
5751 DO i = p_bounds(1, count1), p_bounds(2, count1)
5752 gp(count1) = real(i, kind=
dp)*grid_sp(count1)
5753 k = rec_eval_grid(iw1, ncol, f_vals, count1, gp, grid_sp, step_size, &
5754 istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
5756 ELSE IF (v_count == ncol .AND. ifunc == 1)
THEN
5758 s1v(1, i) = real(i, kind=
dp)*step_size*exp(-lambda*dot_product(gp(:) - f_vals(:, i), &
5759 gp(:) - f_vals(:, i)))
5760 s1v(2, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5765 WRITE (iw1,
'(5F10.5)') gp(:), s1(1)/s1(2)/real(nconf - 1,
dp)
5766 ELSE IF (v_count == ncol .AND. ifunc == 2)
THEN
5768 s1v(1, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5772 WRITE (iw1,
'(5F10.5)') gp(:), -lambda*log(s1(1))
5774 END FUNCTION rec_eval_grid
5787 SUBROUTINE read_frames(frame_section, para_env, nr_frames, r_ref, n_atoms)
5791 INTEGER,
INTENT(IN) :: nr_frames
5792 REAL(
dp),
DIMENSION(:, :),
POINTER :: r_ref
5793 INTEGER,
INTENT(OUT) :: n_atoms
5795 CHARACTER(LEN=default_path_length) :: filename
5796 CHARACTER(LEN=default_string_length) :: dummy_char
5797 INTEGER :: i, j, natom
5798 LOGICAL :: explicit, my_end
5799 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rptr
5812 ALLOCATE (r_ref(3*natom, nr_frames))
5815 cpassert(3*natom ==
SIZE(r_ref, 1))
5819 i_rep_val=j, r_vals=rptr)
5820 r_ref((j - 1)*3 + 1:(j - 1)*3 + 3, i) = rptr(1:3)
5826 cpassert(trim(filename) /=
"")
5828 CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
5834 ALLOCATE (r_ref(3*natom, nr_frames))
5837 cpassert(3*natom ==
SIZE(r_ref, 1))
5843 CALL cp_abort(__location__, &
5844 "Number of lines in XYZ format not equal to the number of atoms."// &
5845 " Error in XYZ format for COORD_A (CV rmsd). Very probably the"// &
5846 " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
5848 READ (parser%input_line, *) dummy_char, rptr(1:3)
5859 END SUBROUTINE read_frames
5870 SUBROUTINE wc_colvar(colvar, cell, subsys, particles, qs_env)
5871 TYPE(colvar_type),
POINTER :: colvar
5872 TYPE(cell_type),
POINTER :: cell
5873 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5874 TYPE(particle_type),
DIMENSION(:), &
5875 OPTIONAL,
POINTER :: particles
5876 TYPE(qs_environment_type),
POINTER,
OPTIONAL :: qs_env
5878 INTEGER :: od, h, oa
5879 REAL(dp) :: rod(3), roa(3), rh(3), &
5880 x, y, s(3), xv(3), dmin, amin
5881 INTEGER :: idmin, iamin, i, j
5882 TYPE(particle_list_type),
POINTER :: particles_i
5883 TYPE(particle_type),
DIMENSION(:), &
5884 POINTER :: my_particles
5885 TYPE(wannier_centres_type),
DIMENSION(:),
POINTER :: wc
5886 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
5887 INTEGER :: nwca, nwcd
5890 NULLIFY (particles_i, wc)
5892 cpassert(colvar%type_id == wc_colvar_id)
5893 IF (
PRESENT(particles))
THEN
5894 my_particles => particles
5896 cpassert(
PRESENT(subsys))
5897 CALL cp_subsys_get(subsys, particles=particles_i)
5898 my_particles => particles_i%els
5900 CALL get_qs_env(qs_env, wanniercentres=wc)
5901 rcut = colvar%Wc%rcut
5902 od = colvar%Wc%ids(1)
5903 h = colvar%Wc%ids(2)
5904 oa = colvar%Wc%ids(3)
5905 CALL get_coordinates(colvar, od, rod, my_particles)
5906 CALL get_coordinates(colvar, h, rh, my_particles)
5907 CALL get_coordinates(colvar, oa, roa, my_particles)
5908 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
5909 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
5912 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
5913 x = distance(rod - wc(1)%centres(:, j))
5914 y = distance(roa - wc(1)%centres(:, j))
5926 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
5927 amin = distance(rh - wc(1)%centres(:, wcai(1)))
5932 x = distance(rh - wc(1)%centres(:, wcdi(i)))
5939 x = distance(rh - wc(1)%centres(:, wcai(i)))
5946 colvar%ss = wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
5956 REAL(dp) function distance(rij)
5957 REAL(dp),
INTENT(in) :: rij(3)
5959 s = matmul(cell%h_inv, rij)
5961 xv = matmul(cell%hmat, s)
5962 distance = norm2(xv)
5963 END FUNCTION distance
5965 END SUBROUTINE wc_colvar
5976 SUBROUTINE hbp_colvar(colvar, cell, subsys, particles, qs_env)
5977 TYPE(colvar_type),
POINTER :: colvar
5978 TYPE(cell_type),
POINTER :: cell
5979 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5980 TYPE(particle_type),
DIMENSION(:), &
5981 OPTIONAL,
POINTER :: particles
5982 TYPE(qs_environment_type),
OPTIONAL,
POINTER :: qs_env
5984 INTEGER :: od, h, oa
5985 REAL(dp) :: rod(3), roa(3), rh(3), &
5986 x, y, s(3), xv(3), dmin, amin
5987 INTEGER :: idmin, iamin, i, j, il, output_unit
5988 TYPE(particle_list_type),
POINTER :: particles_i
5989 TYPE(particle_type),
DIMENSION(:), &
5990 POINTER :: my_particles
5991 TYPE(wannier_centres_type), &
5992 DIMENSION(:),
POINTER :: wc
5993 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
5994 INTEGER :: nwca, nwcd
5997 NULLIFY (particles_i, wc)
5998 output_unit = cp_logger_get_default_io_unit()
6000 cpassert(colvar%type_id == hbp_colvar_id)
6001 IF (
PRESENT(particles))
THEN
6002 my_particles => particles
6004 cpassert(
PRESENT(subsys))
6005 CALL cp_subsys_get(subsys, particles=particles_i)
6006 my_particles => particles_i%els
6008 CALL get_qs_env(qs_env, wanniercentres=wc)
6009 rcut = colvar%HBP%rcut
6010 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
6011 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
6013 DO il = 1, colvar%HBP%nPoints
6014 od = colvar%HBP%ids(il, 1)
6015 h = colvar%HBP%ids(il, 2)
6016 oa = colvar%HBP%ids(il, 3)
6017 CALL get_coordinates(colvar, od, rod, my_particles)
6018 CALL get_coordinates(colvar, h, rh, my_particles)
6019 CALL get_coordinates(colvar, oa, roa, my_particles)
6022 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
6023 x = distance(rod - wc(1)%centres(:, j))
6024 y = distance(roa - wc(1)%centres(:, j))
6036 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
6037 amin = distance(rh - wc(1)%centres(:, wcai(1)))
6042 x = distance(rh - wc(1)%centres(:, wcdi(i)))
6049 x = distance(rh - wc(1)%centres(:, wcai(i)))
6055 colvar%HBP%ewc(il) = colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6056 colvar%ss = colvar%ss + colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6058 IF (output_unit > 0)
THEN
6059 DO il = 1, colvar%HBP%nPoints
6060 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP| = ", colvar%HBP%ewc(il)
6062 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP|\theta(x) = ", colvar%ss
6073 REAL(dp) function distance(rij)
6074 REAL(dp),
INTENT(in) :: rij(3)
6076 s = matmul(cell%h_inv, rij)
6078 xv = matmul(cell%hmat, s)
6079 distance = norm2(xv)
6080 END FUNCTION distance
6082 END SUBROUTINE hbp_colvar
Handles all functions related to the CELL.
subroutine, public cell_transform_input_cartesian(cell, vector)
Transform a Cartesian real-space vector from the user input cell frame into CP2K's canonical internal...
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, cell)
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, cell_ref, use_ref_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, ipi_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...
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
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)
...
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, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
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
Type defining parameters related to the simulation cell.
parameters for a collective variable
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represent a pointer to a subsys, to be able to create arrays of pointers
represents a system: atoms, molecules, their pos,vel,...
wrapper to abstract the force evaluation of the various methods
stores all the informations relevant to an mpi environment
represents a pointer to a list
represent a list of objects