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 :: use_mixed_energy
146 LOGICAL,
DIMENSION(27) :: 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 voronoiipz_section, wc_section
164 TYPE(
section_vals_type),
POINTER :: wrk_section, xyz_diag_section, xyz_outerdiag_section
166 CALL timeset(routinen, handle)
167 NULLIFY (logger, c_kinds, iatms)
169 my_subsection = .false.
177 plane_plane_angle_section &
186 acid_hyd_shell_section &
189 can_return_null=.true.)
190 distance_from_path_section &
192 i_rep_section=icol, can_return_null=.true.)
194 can_return_null=.true.)
203 ring_puckering_section &
217 IF (
ASSOCIATED(reaction_path_section))
THEN
219 explicit=my_subsection(10))
221 IF (
ASSOCIATED(distance_from_path_section))
THEN
223 explicit=my_subsection(16))
225 IF (
ASSOCIATED(combine_section))
THEN
230 explicit=my_subsection(13))
239 explicit=my_subsection(22))
247 cpassert(count(my_subsection) == 1)
248 cpassert(.NOT.
ASSOCIATED(colvar))
250 IF (my_subsection(1))
THEN
252 wrk_section => distance_section
254 CALL colvar_check_points(colvar, distance_section, cell)
256 colvar%dist_param%i_at = iatms(1)
257 colvar%dist_param%j_at = iatms(2)
260 ELSE IF (my_subsection(2))
THEN
262 wrk_section => angle_section
264 CALL colvar_check_points(colvar, angle_section, cell)
266 colvar%angle_param%i_at_angle = iatms
267 ELSE IF (my_subsection(3))
THEN
269 wrk_section => torsion_section
271 CALL colvar_check_points(colvar, torsion_section, cell)
273 colvar%torsion_param%i_at_tors = iatms
274 colvar%torsion_param%o0 = 0.0_dp
275 ELSE IF (my_subsection(4))
THEN
277 wrk_section => coordination_section
279 CALL colvar_check_points(colvar, coordination_section, cell)
280 NULLIFY (colvar%coord_param%i_at_from, colvar%coord_param%c_kinds_from)
281 NULLIFY (colvar%coord_param%i_at_to, colvar%coord_param%c_kinds_to)
282 NULLIFY (colvar%coord_param%i_at_to_b, colvar%coord_param%c_kinds_to_b)
290 CALL reallocate(colvar%coord_param%i_at_from, 1, ndim +
SIZE(iatms))
291 colvar%coord_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
292 ndim = ndim +
SIZE(iatms)
294 colvar%coord_param%n_atoms_from = ndim
295 colvar%coord_param%use_kinds_from = .false.
302 CALL reallocate(colvar%coord_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
303 colvar%coord_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
304 ndim = ndim +
SIZE(c_kinds)
306 colvar%coord_param%n_atoms_from = 0
307 colvar%coord_param%use_kinds_from = .true.
310 CALL uppercase(colvar%coord_param%c_kinds_from(k))
320 CALL reallocate(colvar%coord_param%i_at_to, 1, ndim +
SIZE(iatms))
321 colvar%coord_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
322 ndim = ndim +
SIZE(iatms)
324 colvar%coord_param%n_atoms_to = ndim
325 colvar%coord_param%use_kinds_to = .false.
332 CALL reallocate(colvar%coord_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
333 colvar%coord_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
334 ndim = ndim +
SIZE(c_kinds)
336 colvar%coord_param%n_atoms_to = 0
337 colvar%coord_param%use_kinds_to = .true.
340 CALL uppercase(colvar%coord_param%c_kinds_to(k))
351 IF (n_var /= 0 .OR. n_var_k /= 0)
THEN
352 colvar%coord_param%do_chain = .true.
357 CALL reallocate(colvar%coord_param%i_at_to_b, 1, ndim +
SIZE(iatms))
358 colvar%coord_param%i_at_to_b(ndim + 1:ndim +
SIZE(iatms)) = iatms
359 ndim = ndim +
SIZE(iatms)
361 colvar%coord_param%n_atoms_to_b = ndim
362 colvar%coord_param%use_kinds_to_b = .false.
366 cpassert(n_var_k > 0)
369 CALL reallocate(colvar%coord_param%c_kinds_to_b, 1, ndim +
SIZE(c_kinds))
370 colvar%coord_param%c_kinds_to_b(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
371 ndim = ndim +
SIZE(c_kinds)
373 colvar%coord_param%n_atoms_to_b = 0
374 colvar%coord_param%use_kinds_to_b = .true.
377 CALL uppercase(colvar%coord_param%c_kinds_to_b(k))
385 colvar%coord_param%do_chain = .false.
386 colvar%coord_param%n_atoms_to_b = 0
387 colvar%coord_param%use_kinds_to_b = .false.
388 NULLIFY (colvar%coord_param%i_at_to_b)
389 NULLIFY (colvar%coord_param%c_kinds_to_b)
390 colvar%coord_param%nncrd_b = 0
391 colvar%coord_param%ndcrd_b = 0
392 colvar%coord_param%r_0_b = 0._dp
395 ELSE IF (my_subsection(5))
THEN
397 wrk_section => plane_dist_section
399 CALL colvar_check_points(colvar, plane_dist_section, cell)
401 cpassert(
SIZE(iatms) == 3)
402 colvar%plane_distance_param%plane = iatms
404 colvar%plane_distance_param%point = iatm
406 ELSE IF (my_subsection(6))
THEN
408 wrk_section => rotation_section
410 CALL colvar_check_points(colvar, rotation_section, cell)
411 CALL section_vals_val_get(rotation_section,
"P1_BOND1", i_val=colvar%rotation_param%i_at1_bond1)
412 CALL section_vals_val_get(rotation_section,
"P2_BOND1", i_val=colvar%rotation_param%i_at2_bond1)
413 CALL section_vals_val_get(rotation_section,
"P1_BOND2", i_val=colvar%rotation_param%i_at1_bond2)
414 CALL section_vals_val_get(rotation_section,
"P2_BOND2", i_val=colvar%rotation_param%i_at2_bond2)
415 ELSE IF (my_subsection(7))
THEN
417 wrk_section => dfunct_section
419 CALL colvar_check_points(colvar, dfunct_section, cell)
421 colvar%dfunct_param%i_at_dfunct = iatms
424 ELSE IF (my_subsection(8))
THEN
426 wrk_section => qparm_section
428 CALL colvar_check_points(colvar, qparm_section, cell)
431 CALL section_vals_val_get(qparm_section,
"INCLUDE_IMAGES", l_val=colvar%qparm_param%include_images)
434 NULLIFY (colvar%qparm_param%i_at_from)
435 NULLIFY (colvar%qparm_param%i_at_to)
440 CALL reallocate(colvar%qparm_param%i_at_from, 1, ndim +
SIZE(iatms))
441 colvar%qparm_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
442 ndim = ndim +
SIZE(iatms)
444 colvar%qparm_param%n_atoms_from = ndim
450 CALL reallocate(colvar%qparm_param%i_at_to, 1, ndim +
SIZE(iatms))
451 colvar%qparm_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
452 ndim = ndim +
SIZE(iatms)
454 colvar%qparm_param%n_atoms_to = ndim
455 ELSE IF (my_subsection(9))
THEN
458 NULLIFY (colvar%hydronium_shell_param%i_oxygens)
459 NULLIFY (colvar%hydronium_shell_param%i_hydrogens)
461 colvar%hydronium_shell_param%n_oxygens, &
462 colvar%hydronium_shell_param%n_hydrogens, &
463 colvar%hydronium_shell_param%i_oxygens, &
464 colvar%hydronium_shell_param%i_hydrogens)
465 ELSE IF (my_subsection(10) .OR. my_subsection(16))
THEN
467 IF (my_subsection(10))
THEN
468 path_section => reaction_path_section
472 ELSE IF (my_subsection(16))
THEN
473 path_section => distance_from_path_section
478 colvar%use_points = .false.
480 CALL section_vals_val_get(path_section,
"DISTANCES_RMSD", l_val=colvar%reaction_path_param%dist_rmsd)
482 IF (colvar%reaction_path_param%dist_rmsd .AND. colvar%reaction_path_param%rmsd)
THEN
483 cpabort(
"CV REACTION PATH: only one between DISTANCES_RMSD and RMSD can be used ")
485 IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd)
THEN
486 NULLIFY (colvar%reaction_path_param%i_rmsd, colvar%reaction_path_param%r_ref)
490 colvar%reaction_path_param%nr_frames = nr_frame
491 CALL read_frames(frame_section, para_env, nr_frame, colvar%reaction_path_param%r_ref, &
492 colvar%reaction_path_param%n_components)
494 IF (colvar%reaction_path_param%subset ==
rmsd_all)
THEN
495 ALLOCATE (colvar%reaction_path_param%i_rmsd(colvar%reaction_path_param%n_components))
496 DO i = 1, colvar%reaction_path_param%n_components
497 colvar%reaction_path_param%i_rmsd(i) = i
499 ELSE IF (colvar%reaction_path_param%subset ==
rmsd_list)
THEN
507 CALL reallocate(colvar%reaction_path_param%i_rmsd, 1, ndim +
SIZE(iatms))
508 colvar%reaction_path_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
509 ndim = ndim +
SIZE(iatms)
511 colvar%reaction_path_param%n_components = ndim
513 cpabort(
"CV REACTION PATH: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
517 CALL section_vals_val_get(path_section,
"ALIGN_FRAMES", l_val=colvar%reaction_path_param%align_frames)
521 ALLOCATE (colvar%reaction_path_param%colvar_p(ncol))
524 NULLIFY (colvar%reaction_path_param%colvar_p(i)%colvar)
525 CALL colvar_read(colvar%reaction_path_param%colvar_p(i)%colvar, i, colvar_subsection, para_env, cell)
528 cpabort(
"CV REACTION PATH: the number of CV to define the path must be >0 ")
530 colvar%reaction_path_param%n_components = ncol
533 CALL section_vals_val_get(path_section,
"STEP_SIZE", r_val=colvar%reaction_path_param%step_size)
534 iend = ceiling(max(range(1), range(2))/colvar%reaction_path_param%step_size)
535 istart = floor(min(range(1), range(2))/colvar%reaction_path_param%step_size)
536 colvar%reaction_path_param%function_bounds(1) = istart
537 colvar%reaction_path_param%function_bounds(2) = iend
538 colvar%reaction_path_param%nr_frames = 2
539 ALLOCATE (colvar%reaction_path_param%f_vals(ncol, istart:iend))
543 IF (
SIZE(colvar%reaction_path_param%colvar_p) /= ncol)
THEN
544 CALL cp_abort(__location__, &
545 "Number of defined COLVAR for REACTION_PATH or DISTANCE_FROM_PATH "// &
546 "is different from the number of functions! It is not possible "// &
547 "to define COLVARs in a REACTION_PATH or DISTANCE_FROM_PATH type "// &
548 "COLVAR and omit the corresponding functions!")
553 CALL compress(path_function, full=.true.)
554 CALL parsef(i, trim(path_function), my_par)
556 my_val = real(j, kind=
dp)*colvar%reaction_path_param%step_size
557 colvar%reaction_path_param%f_vals(i, j) =
evalf(i, my_val)
563 "MAP", middle_name=fmid, extension=
".dat", file_status=
"REPLACE")
566 ALLOCATE (grid_sp(ncol))
571 cpassert(ncol ==
SIZE(grid_sp))
572 ALLOCATE (p_range(2, ncol))
573 ALLOCATE (p_bounds(2, ncol))
576 p_range(:, i) = g_range(:)
577 p_bounds(2, i) = ceiling(max(p_range(1, i), p_range(2, i))/grid_sp(i))
578 p_bounds(1, i) = floor(min(p_range(1, i), p_range(2, i))/grid_sp(i))
580 ALLOCATE (s1v(2, istart:iend))
582 ALLOCATE (grid_point(ncol))
584 kk = rec_eval_grid(iw1, ncol, colvar%reaction_path_param%f_vals, v_count, &
585 grid_point, grid_sp, colvar%reaction_path_param%step_size, istart, &
586 iend, s1v, s1, p_bounds, colvar%reaction_path_param%lambda, ifunc=ifunc, &
587 nconf=colvar%reaction_path_param%nr_frames)
590 DEALLOCATE (p_bounds)
593 DEALLOCATE (grid_point)
599 ELSE IF (my_subsection(11))
THEN
602 colvar%use_points = .false.
605 ALLOCATE (colvar%combine_cvs_param%colvar_p(ncol))
608 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
610 WRITE (iw,
'( A )')
' '// &
611 '**********************************************************************'
612 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
613 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| COMBINATION OF THE FOLLOWING COLVARS:'
616 "PRINT%PROGRAM_RUN_INFO")
619 NULLIFY (colvar%combine_cvs_param%colvar_p(i)%colvar)
620 CALL colvar_read(colvar%combine_cvs_param%colvar_p(i)%colvar, i, colvar_subsection, para_env, cell)
624 CALL compress(colvar%combine_cvs_param%function, full=.true.)
627 ALLOCATE (colvar%combine_cvs_param%variables(
SIZE(my_par)))
628 colvar%combine_cvs_param%variables = my_par
630 IF (
SIZE(my_par) /= ncol)
THEN
631 CALL cp_abort(__location__, &
632 "Number of defined COLVAR for COMBINE_COLVAR is different from the "// &
633 "number of variables! It is not possible to define COLVARs in a COMBINE_COLVAR "// &
634 "and avoid their usage in the combining function!")
637 ALLOCATE (colvar%combine_cvs_param%c_parameters(0))
640 isize =
SIZE(colvar%combine_cvs_param%c_parameters)
642 CALL reallocate(colvar%combine_cvs_param%c_parameters, 1, isize +
SIZE(my_par))
643 colvar%combine_cvs_param%c_parameters(isize + 1:isize +
SIZE(my_par)) = my_par
645 ALLOCATE (colvar%combine_cvs_param%v_parameters(0))
648 isize =
SIZE(colvar%combine_cvs_param%v_parameters)
650 CALL reallocate(colvar%combine_cvs_param%v_parameters, 1, isize +
SIZE(my_vals))
651 colvar%combine_cvs_param%v_parameters(isize + 1:isize +
SIZE(my_vals)) = my_vals
656 ELSE IF (my_subsection(12))
THEN
658 wrk_section => population_section
660 CALL colvar_check_points(colvar, population_section, cell)
662 NULLIFY (colvar%population_param%i_at_from, colvar%population_param%c_kinds_from)
663 NULLIFY (colvar%population_param%i_at_to, colvar%population_param%c_kinds_to)
672 CALL reallocate(colvar%population_param%i_at_from, 1, ndim +
SIZE(iatms))
673 colvar%population_param%i_at_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
674 ndim = ndim +
SIZE(iatms)
676 colvar%population_param%n_atoms_from = ndim
677 colvar%population_param%use_kinds_from = .false.
684 CALL reallocate(colvar%population_param%c_kinds_from, 1, ndim +
SIZE(c_kinds))
685 colvar%population_param%c_kinds_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
686 ndim = ndim +
SIZE(c_kinds)
688 colvar%population_param%n_atoms_from = 0
689 colvar%population_param%use_kinds_from = .true.
692 CALL uppercase(colvar%population_param%c_kinds_from(k))
702 CALL reallocate(colvar%population_param%i_at_to, 1, ndim +
SIZE(iatms))
703 colvar%population_param%i_at_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
704 ndim = ndim +
SIZE(iatms)
706 colvar%population_param%n_atoms_to = ndim
707 colvar%population_param%use_kinds_to = .false.
714 CALL reallocate(colvar%population_param%c_kinds_to, 1, ndim +
SIZE(c_kinds))
715 colvar%population_param%c_kinds_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
716 ndim = ndim +
SIZE(c_kinds)
718 colvar%population_param%n_atoms_to = 0
719 colvar%population_param%use_kinds_to = .true.
722 CALL uppercase(colvar%population_param%c_kinds_to(k))
731 ELSE IF (my_subsection(13))
THEN
733 wrk_section => plane_plane_angle_section
735 CALL colvar_check_points(colvar, plane_plane_angle_section, cell)
740 cpabort(
"PLANE_PLANE_ANGLE Colvar section: Two PLANE sections must be provided!")
744 i_val=colvar%plane_plane_angle_param%plane1%type_of_def)
745 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_vec)
THEN
748 colvar%plane_plane_angle_param%plane1%normal_vec = s1
749 IF (
PRESENT(cell))
THEN
750 IF (
ASSOCIATED(cell))
THEN
757 colvar%plane_plane_angle_param%plane1%points = iatms
762 i_val=colvar%plane_plane_angle_param%plane2%type_of_def)
763 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_vec)
THEN
766 colvar%plane_plane_angle_param%plane2%normal_vec = s1
767 IF (
PRESENT(cell))
THEN
768 IF (
ASSOCIATED(cell))
THEN
775 colvar%plane_plane_angle_param%plane2%points = iatms
777 ELSE IF (my_subsection(14))
THEN
779 wrk_section => gyration_section
781 CALL colvar_check_points(colvar, gyration_section, cell)
783 NULLIFY (colvar%gyration_param%i_at, colvar%gyration_param%c_kinds)
792 CALL reallocate(colvar%gyration_param%i_at, 1, ndim +
SIZE(iatms))
793 colvar%gyration_param%i_at(ndim + 1:ndim +
SIZE(iatms)) = iatms
794 ndim = ndim +
SIZE(iatms)
796 colvar%gyration_param%n_atoms = ndim
797 colvar%gyration_param%use_kinds = .false.
804 CALL reallocate(colvar%gyration_param%c_kinds, 1, ndim +
SIZE(c_kinds))
805 colvar%gyration_param%c_kinds(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
806 ndim = ndim +
SIZE(c_kinds)
808 colvar%gyration_param%n_atoms = 0
809 colvar%gyration_param%use_kinds = .true.
812 CALL uppercase(colvar%gyration_param%c_kinds(k))
815 ELSE IF (my_subsection(15))
THEN
817 wrk_section => rmsd_section
820 NULLIFY (colvar%rmsd_param%i_rmsd, colvar%rmsd_param%r_ref, colvar%rmsd_param%weights)
825 colvar%rmsd_param%nr_frames = nr_frame
827 cpassert(nr_frame >= 1 .AND. nr_frame <= 2)
828 CALL read_frames(frame_section, para_env, nr_frame, colvar%rmsd_param%r_ref, &
829 colvar%rmsd_param%n_atoms)
830 ALLOCATE (colvar%rmsd_param%weights(colvar%rmsd_param%n_atoms))
831 colvar%rmsd_param%weights = 0.0_dp
833 IF (colvar%rmsd_param%subset ==
rmsd_all)
THEN
834 ALLOCATE (colvar%rmsd_param%i_rmsd(colvar%rmsd_param%n_atoms))
835 DO i = 1, colvar%rmsd_param%n_atoms
836 colvar%rmsd_param%i_rmsd(i) = i
838 ELSE IF (colvar%rmsd_param%subset ==
rmsd_list)
THEN
846 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
847 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
848 ndim = ndim +
SIZE(iatms)
850 colvar%rmsd_param%n_atoms = ndim
852 cpabort(
"CV RMSD: if SUBSET_TYPE=LIST a list of atoms needs to be provided ")
861 CALL reallocate(colvar%rmsd_param%i_rmsd, 1, ndim +
SIZE(iatms))
862 colvar%rmsd_param%i_rmsd(ndim + 1:ndim +
SIZE(iatms)) = iatms
863 ndim = ndim +
SIZE(iatms)
865 colvar%rmsd_param%n_atoms = ndim
867 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of atoms needs to be provided ")
876 weights(ndim + 1:ndim +
SIZE(wei)) = wei
877 ndim = ndim +
SIZE(wei)
879 IF (ndim /= colvar%rmsd_param%n_atoms)
THEN
880 CALL cp_abort(__location__,
"CV RMSD: list of atoms and list of "// &
881 "weights need to contain same number of entries. ")
884 ii = colvar%rmsd_param%i_rmsd(i)
885 colvar%rmsd_param%weights(ii) = weights(i)
889 cpabort(
"CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of weights need to be provided. ")
893 cpabort(
"CV RMSD: unknown SUBSET_TYPE.")
898 ELSE IF (my_subsection(17))
THEN
900 wrk_section => xyz_diag_section
902 CALL colvar_check_points(colvar, wrk_section, cell)
906 CALL section_vals_val_get(wrk_section,
"ABSOLUTE_POSITION", l_val=colvar%xyz_diag_param%use_absolute_position)
907 colvar%xyz_diag_param%i_atom = iatm
908 colvar%xyz_diag_param%component = icomponent
909 ELSE IF (my_subsection(18))
THEN
911 wrk_section => xyz_outerdiag_section
913 CALL colvar_check_points(colvar, wrk_section, cell)
915 colvar%xyz_outerdiag_param%i_atoms = iatms
917 colvar%xyz_outerdiag_param%components(1) = icomponent
919 colvar%xyz_outerdiag_param%components(2) = icomponent
921 ELSE IF (my_subsection(19))
THEN
923 wrk_section => u_section
926 CALL section_vals_get(colvar%u_param%mixed_energy_section, explicit=use_mixed_energy)
927 IF (.NOT. use_mixed_energy)
NULLIFY (colvar%u_param%mixed_energy_section)
928 ELSE IF (my_subsection(20))
THEN
930 wrk_section => wc_section
932 CALL colvar_check_points(colvar, wc_section, cell)
936 colvar%Wc%ids = iatms
937 ELSE IF (my_subsection(21))
THEN
939 wrk_section => hbp_section
941 CALL colvar_check_points(colvar, hbp_section, cell)
947 ALLOCATE (colvar%HBP%ids(colvar%HBP%nPoints, 3))
948 ALLOCATE (colvar%HBP%ewc(colvar%HBP%nPoints))
949 DO i = 1, colvar%HBP%nPoints
951 colvar%HBP%ids(i, :) = iatms
953 ELSE IF (my_subsection(22))
THEN
957 colvar%ring_puckering_param%nring =
SIZE(iatms)
958 ALLOCATE (colvar%ring_puckering_param%atoms(
SIZE(iatms)))
959 colvar%ring_puckering_param%atoms = iatms
961 i_val=colvar%ring_puckering_param%iq)
963 ndim = colvar%ring_puckering_param%nring
965 cpabort(
"CV Ring Puckering: Ring size has to be 4 or larger. ")
967 ii = colvar%ring_puckering_param%iq
968 IF (abs(ii) == 1 .OR. ii < -(ndim - 1)/2 .OR. ii > ndim/2)
THEN
969 cpabort(
"CV Ring Puckering: Invalid coordinate number.")
971 ELSE IF (my_subsection(23))
THEN
973 wrk_section => mindist_section
975 CALL colvar_check_points(colvar, mindist_section, cell)
976 NULLIFY (colvar%mindist_param%i_dist_from, colvar%mindist_param%i_coord_from, &
977 colvar%mindist_param%k_coord_from, colvar%mindist_param%i_coord_to, &
978 colvar%mindist_param%k_coord_to)
980 colvar%mindist_param%n_dist_from =
SIZE(iatms)
981 ALLOCATE (colvar%mindist_param%i_dist_from(
SIZE(iatms)))
982 colvar%mindist_param%i_dist_from = iatms
989 CALL reallocate(colvar%mindist_param%i_coord_from, 1, ndim +
SIZE(iatms))
990 colvar%mindist_param%i_coord_from(ndim + 1:ndim +
SIZE(iatms)) = iatms
991 ndim = ndim +
SIZE(iatms)
993 colvar%mindist_param%n_coord_from = ndim
994 colvar%mindist_param%use_kinds_from = .false.
1001 CALL reallocate(colvar%mindist_param%k_coord_from, 1, ndim +
SIZE(c_kinds))
1002 colvar%mindist_param%k_coord_from(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
1003 ndim = ndim +
SIZE(c_kinds)
1005 colvar%mindist_param%n_coord_from = 0
1006 colvar%mindist_param%use_kinds_from = .true.
1009 CALL uppercase(colvar%mindist_param%k_coord_from(k))
1015 IF (n_var /= 0)
THEN
1019 CALL reallocate(colvar%mindist_param%i_coord_to, 1, ndim +
SIZE(iatms))
1020 colvar%mindist_param%i_coord_to(ndim + 1:ndim +
SIZE(iatms)) = iatms
1021 ndim = ndim +
SIZE(iatms)
1023 colvar%mindist_param%n_coord_to = ndim
1024 colvar%mindist_param%use_kinds_to = .false.
1031 CALL reallocate(colvar%mindist_param%k_coord_to, 1, ndim +
SIZE(c_kinds))
1032 colvar%mindist_param%k_coord_to(ndim + 1:ndim +
SIZE(c_kinds)) = c_kinds
1033 ndim = ndim +
SIZE(c_kinds)
1035 colvar%mindist_param%n_coord_to = 0
1036 colvar%mindist_param%use_kinds_to = .true.
1039 CALL uppercase(colvar%mindist_param%k_coord_to(k))
1048 ELSE IF (my_subsection(24))
THEN
1051 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_water)
1052 NULLIFY (colvar%acid_hyd_dist_param%i_oxygens_acid)
1053 NULLIFY (colvar%acid_hyd_dist_param%i_hydrogens)
1055 colvar%acid_hyd_dist_param%n_oxygens_water, &
1056 colvar%acid_hyd_dist_param%n_oxygens_acid, &
1057 colvar%acid_hyd_dist_param%n_hydrogens, &
1058 colvar%acid_hyd_dist_param%i_oxygens_water, &
1059 colvar%acid_hyd_dist_param%i_oxygens_acid, &
1060 colvar%acid_hyd_dist_param%i_hydrogens)
1061 ELSE IF (my_subsection(25))
THEN
1064 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_water)
1065 NULLIFY (colvar%acid_hyd_shell_param%i_oxygens_acid)
1066 NULLIFY (colvar%acid_hyd_shell_param%i_hydrogens)
1068 colvar%acid_hyd_shell_param%n_oxygens_water, &
1069 colvar%acid_hyd_shell_param%n_oxygens_acid, &
1070 colvar%acid_hyd_shell_param%n_hydrogens, &
1071 colvar%acid_hyd_shell_param%i_oxygens_water, &
1072 colvar%acid_hyd_shell_param%i_oxygens_acid, &
1073 colvar%acid_hyd_shell_param%i_hydrogens)
1074 ELSE IF (my_subsection(26))
THEN
1077 NULLIFY (colvar%hydronium_dist_param%i_oxygens)
1078 NULLIFY (colvar%hydronium_dist_param%i_hydrogens)
1080 colvar%hydronium_dist_param%n_oxygens, &
1081 colvar%hydronium_dist_param%n_hydrogens, &
1082 colvar%hydronium_dist_param%i_oxygens, &
1083 colvar%hydronium_dist_param%i_hydrogens)
1084 ELSE IF (my_subsection(27))
THEN
1086 cpassert(
ASSOCIATED(voronoiipz_section))
1087 CALL read_voronoiipz_colvars(voronoiipz_section, colvar)
1092 "PRINT%PROGRAM_RUN_INFO", extension=
".colvarLog")
1095 IF (colvar%use_points) tag =
"POINTS:"
1098 WRITE (iw,
'( A )')
' '// &
1099 '----------------------------------------------------------------------'
1100 WRITE (iw,
'( A,I8)')
' COLVARS| COLVAR INPUT INDEX: ', icol
1103 SELECT CASE (colvar%type_id)
1105 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| ANGLE >>> '//tag, &
1106 colvar%angle_param%i_at_angle
1108 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| DISTANCE DIFFERENCE >>> '//tag, &
1109 colvar%dfunct_param%i_at_dfunct
1111 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE DISTANCE - PLANE >>> '//tag, &
1112 colvar%plane_distance_param%plane
1113 WRITE (iw,
'( A,T73,1I8)')
' COLVARS| PLANE DISTANCE - POINT >>> '//tag, &
1114 colvar%plane_distance_param%point
1116 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
1117 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1118 colvar%plane_plane_angle_param%plane1%points
1120 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1121 colvar%plane_plane_angle_param%plane1%normal_vec
1124 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
1125 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag, &
1126 colvar%plane_plane_angle_param%plane2%points
1128 WRITE (iw,
'( A,T57,3F8.3)')
' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag, &
1129 colvar%plane_plane_angle_param%plane2%normal_vec
1132 WRITE (iw,
'( A,T49,4I8)')
' COLVARS| TORSION >>> '//tag, &
1133 colvar%torsion_param%i_at_tors
1135 WRITE (iw,
'( A,T65,2I8)')
' COLVARS| BOND >>> '//tag, &
1136 colvar%dist_param%i_at, colvar%dist_param%j_at
1138 IF (colvar%coord_param%do_chain)
THEN
1139 WRITE (iw,
'( A)')
' COLVARS| COORDINATION CHAIN FC(from->to)*FC(to->to_B)>> '
1141 IF (colvar%coord_param%use_kinds_from)
THEN
1142 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> FROM KINDS', &
1143 adjustr(colvar%coord_param%c_kinds_from(kk) (1:10)), &
1144 kk=1,
SIZE(colvar%coord_param%c_kinds_from))
1146 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> FROM '//tag, &
1147 colvar%coord_param%i_at_from(kk), &
1148 kk=1,
SIZE(colvar%coord_param%i_at_from))
1150 IF (colvar%coord_param%use_kinds_to)
THEN
1151 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS', &
1152 adjustr(colvar%coord_param%c_kinds_to(kk) (1:10)), &
1153 kk=1,
SIZE(colvar%coord_param%c_kinds_to))
1155 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag, &
1156 colvar%coord_param%i_at_to(kk), &
1157 kk=1,
SIZE(colvar%coord_param%i_at_to))
1159 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%coord_param%r_0
1160 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%coord_param%nncrd
1161 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%coord_param%ndcrd
1162 IF (colvar%coord_param%do_chain)
THEN
1163 IF (colvar%coord_param%use_kinds_to_b)
THEN
1164 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COORDINATION >>> TO KINDS B', &
1165 adjustr(colvar%coord_param%c_kinds_to_b(kk) (1:10)), &
1166 kk=1,
SIZE(colvar%coord_param%c_kinds_to_b))
1168 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COORDINATION >>> TO '//tag//
' B', &
1169 colvar%coord_param%i_at_to_b(kk), &
1170 kk=1,
SIZE(colvar%coord_param%i_at_to_b))
1172 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0 B', colvar%coord_param%r_0_b
1173 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN B', colvar%coord_param%nncrd_b
1174 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND B', colvar%coord_param%ndcrd_b
1177 IF (colvar%population_param%use_kinds_from)
THEN
1178 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> FROM KINDS', &
1179 adjustr(colvar%population_param%c_kinds_from(kk) (1:10)), &
1180 kk=1,
SIZE(colvar%population_param%c_kinds_from))
1182 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> FROM '//tag, &
1183 colvar%population_param%i_at_from(kk), &
1184 kk=1,
SIZE(colvar%population_param%i_at_from))
1186 IF (colvar%population_param%use_kinds_to)
THEN
1187 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| POPULATION based on coordination >>> TO KINDS', &
1188 adjustr(colvar%population_param%c_kinds_to(kk) (1:10)), &
1189 kk=1,
SIZE(colvar%population_param%c_kinds_to))
1191 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| POPULATION based on coordination >>> TO '//tag, &
1192 colvar%population_param%i_at_to(kk), &
1193 kk=1,
SIZE(colvar%population_param%i_at_to))
1195 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%population_param%r_0
1196 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%population_param%nncrd
1197 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%population_param%ndcrd
1198 WRITE (iw,
'( A,T71,I10)')
' COLVARS| N0', colvar%population_param%n0
1199 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| SIGMA', colvar%population_param%sigma
1201 IF (colvar%gyration_param%use_kinds)
THEN
1202 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| Gyration Radius >>> KINDS', &
1203 adjustr(colvar%gyration_param%c_kinds(kk) (1:10)), &
1204 kk=1,
SIZE(colvar%gyration_param%c_kinds))
1206 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Gyration Radius >>> ATOMS '//tag, &
1207 colvar%gyration_param%i_at(kk), &
1208 kk=1,
SIZE(colvar%gyration_param%i_at))
1211 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 1 >>> '//tag, &
1212 colvar%rotation_param%i_at1_bond1
1213 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 1 >>> '//tag, &
1214 colvar%rotation_param%i_at2_bond1
1215 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 1 LINE 2 >>> '//tag, &
1216 colvar%rotation_param%i_at1_bond2
1217 WRITE (iw,
'( A,T71,I10)')
' COLVARS| BOND_ROTATION - POINT 2 LINE 2 >>> '//tag, &
1218 colvar%rotation_param%i_at2_bond2
1220 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> FROM '//tag, &
1221 colvar%qparm_param%i_at_from(kk), &
1222 kk=1,
SIZE(colvar%qparm_param%i_at_from))
1223 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| Q-PARM >>> TO '//tag, &
1224 colvar%qparm_param%i_at_to(kk), &
1225 kk=1,
SIZE(colvar%qparm_param%i_at_to))
1226 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RCUT', colvar%qparm_param%rcut
1227 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RSTART', colvar%qparm_param%rstart
1228 WRITE (iw,
'( A,T71,L10)')
' COLVARS| INCLUDE IMAGES', colvar%qparm_param%include_images
1230 WRITE (iw,
'( A,T71,I10)')
' COLVARS| L', colvar%qparm_param%l
1232 WRITE (iw,
'( A)')
' COLVARS| COMBINING FUNCTION : '// &
1233 trim(colvar%combine_cvs_param%function)
1234 WRITE (iw,
'( A)', advance=
"NO")
' COLVARS| VARIABLES : '
1235 DO i = 1,
SIZE(colvar%combine_cvs_param%variables)
1236 WRITE (iw,
'( A)', advance=
"NO") &
1237 trim(colvar%combine_cvs_param%variables(i))//
" "
1240 WRITE (iw,
'( A)')
' COLVARS| DEFINED PARAMETERS [label] [value]:'
1241 DO i = 1,
SIZE(colvar%combine_cvs_param%c_parameters)
1242 WRITE (iw,
'( A,A7,F9.3)')
' ', &
1243 trim(colvar%combine_cvs_param%c_parameters(i)), colvar%combine_cvs_param%v_parameters(i)
1245 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| ERROR ON DERIVATIVE EVALUATION', &
1246 colvar%combine_cvs_param%lerr
1247 WRITE (iw,
'( A,T71,G10.5)')
' COLVARS| DX', &
1248 colvar%combine_cvs_param%dx
1250 cpwarn(
"Description header for REACTION_PATH COLVAR missing!")
1252 cpwarn(
"Description header for REACTION_PATH COLVAR missing!")
1254 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_shell_param%poh
1255 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_shell_param%qoh
1256 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%hydronium_shell_param%poo
1257 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%hydronium_shell_param%qoo
1258 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%hydronium_shell_param%roo
1259 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_shell_param%roh
1260 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_shell_param%nh
1261 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%hydronium_shell_param%lambda
1263 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POH', colvar%hydronium_dist_param%poh
1264 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOH', colvar%hydronium_dist_param%qoh
1265 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROH', colvar%hydronium_dist_param%roh
1266 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%hydronium_dist_param%pm
1267 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%hydronium_dist_param%qm
1268 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%hydronium_dist_param%nh
1269 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PF', colvar%hydronium_dist_param%pf
1270 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QF', colvar%hydronium_dist_param%qf
1271 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NN', colvar%hydronium_dist_param%nn
1273 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_dist_param%paoh
1274 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_dist_param%qaoh
1275 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_dist_param%pwoh
1276 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_dist_param%qwoh
1277 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_dist_param%pcut
1278 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_dist_param%qcut
1279 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_dist_param%raoh
1280 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_dist_param%rwoh
1281 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_dist_param%nc
1282 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_dist_param%lambda
1284 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PAOH', colvar%acid_hyd_shell_param%paoh
1285 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QAOH', colvar%acid_hyd_shell_param%qaoh
1286 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PWOH', colvar%acid_hyd_shell_param%pwoh
1287 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QWOH', colvar%acid_hyd_shell_param%qwoh
1288 WRITE (iw,
'( A,T71,I10)')
' COLVARS| POO', colvar%acid_hyd_shell_param%poo
1289 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QOO', colvar%acid_hyd_shell_param%qoo
1290 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PM', colvar%acid_hyd_shell_param%pm
1291 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QM', colvar%acid_hyd_shell_param%qm
1292 WRITE (iw,
'( A,T71,I10)')
' COLVARS| PCUT', colvar%acid_hyd_shell_param%pcut
1293 WRITE (iw,
'( A,T71,I10)')
' COLVARS| QCUT', colvar%acid_hyd_shell_param%qcut
1294 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RAOH', colvar%acid_hyd_shell_param%raoh
1295 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| RWOH', colvar%acid_hyd_shell_param%rwoh
1296 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| ROO', colvar%acid_hyd_shell_param%roo
1297 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NH', colvar%acid_hyd_shell_param%nh
1298 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| NC', colvar%acid_hyd_shell_param%nc
1299 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%acid_hyd_shell_param%lambda
1301 cpwarn(
"Description header for RMSD COLVAR missing!")
1303 NULLIFY (section, keyword, enum)
1307 tag_comp = trim(
enum_i2c(enum, colvar%xyz_diag_param%component))
1310 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| POSITION ('//trim(tag_comp) &
1311 //
') >>> '//tag, colvar%xyz_diag_param%i_atom
1313 NULLIFY (section, keyword, enum)
1317 tag_comp1 = trim(
enum_i2c(enum, colvar%xyz_outerdiag_param%components(1)))
1320 tag_comp2 = trim(
enum_i2c(enum, colvar%xyz_outerdiag_param%components(2)))
1323 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| CROSS TERM POSITION ('//trim(tag_comp1) &
1324 //
" * "//trim(tag_comp2)//
') >>> '//tag, colvar%xyz_outerdiag_param%i_atoms
1326 WRITE (iw,
'( A,T77,A4)')
' COLVARS| ENERGY >>> '//tag,
'all!'
1328 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| Wc >>> RCUT: ', &
1330 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| Wc >>> '//tag, &
1333 WRITE (iw,
'( A,T57,I8)')
' COLVARS| HBP >>> NPOINTS', &
1335 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1337 WRITE (iw,
'( A,T57,F16.8)')
' COLVARS| HBP >>> RCUT', &
1339 DO i = 1, colvar%HBP%nPoints
1340 WRITE (iw,
'( A,T57,3I8)')
' COLVARS| HBP >>> '//tag, &
1341 colvar%HBP%ids(i, :)
1344 WRITE (iw,
'( A,T57,I8)')
' COLVARS| Ring Puckering >>> ring size', &
1345 colvar%ring_puckering_param%nring
1346 IF (colvar%ring_puckering_param%iq == 0)
THEN
1347 WRITE (iw,
'( A,T40,A)')
' COLVARS| Ring Puckering >>> coordinate', &
1348 ' Total Puckering Amplitude'
1349 ELSE IF (colvar%ring_puckering_param%iq > 0)
THEN
1350 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1351 ' Puckering Amplitude', &
1352 colvar%ring_puckering_param%iq
1354 WRITE (iw,
'( A,T35,A,T57,I8)')
' COLVARS| Ring Puckering >>> coordinate', &
1355 ' Puckering Angle', &
1356 colvar%ring_puckering_param%iq
1359 WRITE (iw,
'( A)')
' COLVARS| CONDITIONED DISTANCE>> '
1360 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DISTANCE >>> DISTANCE FROM '//tag, &
1361 colvar%mindist_param%i_dist_from(kk), &
1362 kk=1,
SIZE(colvar%mindist_param%i_dist_from))
1363 IF (colvar%mindist_param%use_kinds_from)
THEN
1364 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM KINDS ', &
1365 adjustr(colvar%mindist_param%k_coord_from(kk) (1:10)), &
1366 kk=1,
SIZE(colvar%mindist_param%k_coord_from))
1368 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION FROM '//tag, &
1369 colvar%mindist_param%i_coord_from(kk), &
1370 kk=1,
SIZE(colvar%mindist_param%i_coord_from))
1372 IF (colvar%mindist_param%use_kinds_to)
THEN
1373 WRITE (iw,
'( A,T71,A10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO KINDS ', &
1374 adjustr(colvar%mindist_param%k_coord_to(kk) (1:10)), &
1375 kk=1,
SIZE(colvar%mindist_param%k_coord_to))
1377 WRITE (iw,
'( A,T71,I10)') (
' COLVARS| COND.DIST. >>> COORDINATION TO '//tag, &
1378 colvar%mindist_param%i_coord_to(kk), &
1379 kk=1,
SIZE(colvar%mindist_param%i_coord_to))
1381 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| R0', colvar%mindist_param%r_cut
1382 WRITE (iw,
'( A,T71,I10)')
' COLVARS| NN', colvar%mindist_param%p_exp
1383 WRITE (iw,
'( A,T71,I10)')
' COLVARS| ND', colvar%mindist_param%q_exp
1384 WRITE (iw,
'( A,T71,F10.5)')
' COLVARS| LAMBDA', colvar%mindist_param%lambda
1387 IF (colvar%use_points)
THEN
1388 WRITE (iw,
'( A)')
' COLVARS| INFORMATION ON DEFINED GEOMETRICAL POINTS'
1389 DO kk = 1,
SIZE(colvar%points)
1393 WRITE (iw,
'( A)')
' COLVARS| POINT Nr.'//trim(tmpstr2)//
' OF TYPE: '//trim(tmpstr)
1394 IF (
ASSOCIATED(colvar%points(kk)%atoms))
THEN
1395 WRITE (iw,
'( A)')
' COLVARS| ATOMS BUILDING THE GEOMETRICAL POINT'
1396 WRITE (iw,
'( A, I10)') (
' COLVARS| ATOM:', colvar%points(kk)%atoms(k), k=1,
SIZE(colvar%points(kk)%atoms))
1398 WRITE (iw,
'( A,4X,3F12.6)')
' COLVARS| XYZ POSITION OF FIXED POINT:', colvar%points(kk)%r
1404 WRITE (iw,
'( A )')
' '// &
1405 '----------------------------------------------------------------------'
1407 WRITE (iw,
'( A )')
' '// &
1408 '**********************************************************************'
1412 "PRINT%PROGRAM_RUN_INFO")
1413 CALL timestop(handle)
1427 SUBROUTINE read_hydronium_colvars(section, colvar, colvar_id, n_oxygens, n_hydrogens, &
1428 i_oxygens, i_hydrogens)
1431 INTEGER,
INTENT(IN) :: colvar_id
1432 INTEGER,
INTENT(OUT) :: n_oxygens, n_hydrogens
1433 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens, i_hydrogens
1435 INTEGER :: k, n_var, ndim
1436 INTEGER,
DIMENSION(:),
POINTER :: iatms
1444 CALL reallocate(i_oxygens, 1, ndim +
SIZE(iatms))
1445 i_oxygens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1446 ndim = ndim +
SIZE(iatms)
1454 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1455 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1456 ndim = ndim +
SIZE(iatms)
1460 SELECT CASE (colvar_id)
1485 END SUBROUTINE read_hydronium_colvars
1501 SUBROUTINE read_acid_hydronium_colvars(section, colvar, colvar_id, n_oxygens_water, &
1502 n_oxygens_acid, n_hydrogens, i_oxygens_water, &
1503 i_oxygens_acid, i_hydrogens)
1506 INTEGER,
INTENT(IN) :: colvar_id
1507 INTEGER,
INTENT(OUT) :: n_oxygens_water, n_oxygens_acid, &
1509 INTEGER,
DIMENSION(:),
POINTER :: i_oxygens_water, i_oxygens_acid, &
1512 INTEGER :: k, n_var, ndim
1513 INTEGER,
DIMENSION(:),
POINTER :: iatms
1521 CALL reallocate(i_oxygens_water, 1, ndim +
SIZE(iatms))
1522 i_oxygens_water(ndim + 1:ndim +
SIZE(iatms)) = iatms
1523 ndim = ndim +
SIZE(iatms)
1525 n_oxygens_water = ndim
1531 CALL reallocate(i_oxygens_acid, 1, ndim +
SIZE(iatms))
1532 i_oxygens_acid(ndim + 1:ndim +
SIZE(iatms)) = iatms
1533 ndim = ndim +
SIZE(iatms)
1535 n_oxygens_acid = ndim
1541 CALL reallocate(i_hydrogens, 1, ndim +
SIZE(iatms))
1542 i_hydrogens(ndim + 1:ndim +
SIZE(iatms)) = iatms
1543 ndim = ndim +
SIZE(iatms)
1547 SELECT CASE (colvar_id)
1578 END SUBROUTINE read_acid_hydronium_colvars
1587 SUBROUTINE colvar_check_points(colvar, section, cell)
1590 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
1592 INTEGER :: i, irep, natoms, npoints, nrep, nweights
1593 INTEGER,
DIMENSION(:),
POINTER :: atoms
1595 REAL(kind=
dp),
DIMENSION(:),
POINTER :: r, weights
1598 NULLIFY (point_sections)
1601 cpassert(
ASSOCIATED(colvar))
1605 colvar%use_points = .true.
1607 ALLOCATE (colvar%points(npoints))
1612 NULLIFY (colvar%points(i)%atoms)
1613 NULLIFY (colvar%points(i)%weights)
1614 CALL section_vals_val_get(point_sections,
"TYPE", i_rep_section=i, i_val=colvar%points(i)%type_id)
1615 SELECT CASE (colvar%points(i)%type_id)
1618 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, n_rep_val=nrep, i_vals=atoms)
1620 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1621 natoms = natoms +
SIZE(atoms)
1623 ALLOCATE (colvar%points(i)%atoms(natoms))
1626 CALL section_vals_val_get(point_sections,
"ATOMS", i_rep_section=i, i_rep_val=irep, i_vals=atoms)
1627 colvar%points(i)%atoms(natoms + 1:) = atoms(:)
1628 natoms = natoms +
SIZE(atoms)
1631 ALLOCATE (colvar%points(i)%weights(natoms))
1632 colvar%points(i)%weights = 1.0_dp/real(natoms, kind=
dp)
1638 colvar%points(i)%weights(nweights + 1:) = weights(:)
1639 nweights = nweights +
SIZE(weights)
1641 cpassert(natoms == nweights)
1646 colvar%points(i)%r = r
1647 IF (
PRESENT(cell))
THEN
1653 END SUBROUTINE colvar_check_points
1669 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
1672 OPTIONAL,
POINTER :: fixd_list
1675 LOGICAL :: colvar_ok
1677 colvar_ok =
ASSOCIATED(colvar)
1680 IF (
PRESENT(pos))
THEN
1681 DO i = 1,
SIZE(colvar%i_atom)
1682 j = colvar%i_atom(i)
1683 particles(j)%r = pos(:, j)
1687 colvar%dsdr = 0.0_dp
1688 SELECT CASE (colvar%type_id)
1690 CALL dist_colvar(colvar, cell, particles=particles)
1692 CALL coord_colvar(colvar, cell, particles=particles)
1694 CALL population_colvar(colvar, cell, particles=particles)
1696 CALL gyration_radius_colvar(colvar, cell, particles=particles)
1698 CALL torsion_colvar(colvar, cell, particles=particles)
1700 CALL angle_colvar(colvar, cell, particles=particles)
1702 CALL dfunct_colvar(colvar, cell, particles=particles)
1704 CALL plane_distance_colvar(colvar, cell, particles=particles)
1706 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1708 CALL rotation_colvar(colvar, cell, particles=particles)
1710 CALL qparm_colvar(colvar, cell, particles=particles)
1712 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1714 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1716 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1718 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1720 CALL rmsd_colvar(colvar, particles=particles)
1722 CALL reaction_path_colvar(colvar, cell, particles=particles)
1724 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1726 CALL combine_colvar(colvar, cell, particles=particles)
1728 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1730 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1732 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1734 CALL mindist_colvar(colvar, cell, particles=particles)
1736 cpabort(
"need force_env!")
1739 CALL wc_colvar(colvar, cell, particles=particles)
1742 CALL hbp_colvar(colvar, cell, particles=particles)
1744 CALL colvar_eval_voronoiipz(colvar, cell, particles=particles)
1746 cpabort(
"Unknown colvar type for colvar_eval_mol_f")
1766 LOGICAL :: colvar_ok
1772 NULLIFY (subsys, cell, colvar, qs_env)
1773 CALL force_env_get(force_env, subsys=subsys, cell=cell, qs_env=qs_env)
1774 colvar_ok =
ASSOCIATED(subsys%colvar_p)
1777 colvar => subsys%colvar_p(icolvar)%colvar
1779 colvar%dsdr = 0.0_dp
1780 SELECT CASE (colvar%type_id)
1782 CALL dist_colvar(colvar, cell, subsys=subsys)
1784 CALL coord_colvar(colvar, cell, subsys=subsys)
1786 CALL population_colvar(colvar, cell, subsys=subsys)
1788 CALL gyration_radius_colvar(colvar, cell, subsys=subsys)
1790 CALL torsion_colvar(colvar, cell, subsys=subsys, no_riemann_sheet_op=.true.)
1792 CALL angle_colvar(colvar, cell, subsys=subsys)
1794 CALL dfunct_colvar(colvar, cell, subsys=subsys)
1796 CALL plane_distance_colvar(colvar, cell, subsys=subsys)
1798 CALL plane_plane_angle_colvar(colvar, cell, subsys=subsys)
1800 CALL rotation_colvar(colvar, cell, subsys=subsys)
1802 CALL qparm_colvar(colvar, cell, subsys=subsys)
1804 CALL hydronium_shell_colvar(colvar, cell, subsys=subsys)
1806 CALL hydronium_dist_colvar(colvar, cell, subsys=subsys)
1808 CALL acid_hyd_dist_colvar(colvar, cell, subsys=subsys)
1810 CALL acid_hyd_shell_colvar(colvar, cell, subsys=subsys)
1812 CALL rmsd_colvar(colvar, subsys=subsys)
1814 CALL reaction_path_colvar(colvar, cell, subsys=subsys)
1816 CALL distance_from_path_colvar(colvar, cell, subsys=subsys)
1818 CALL combine_colvar(colvar, cell, subsys=subsys)
1820 CALL xyz_diag_colvar(colvar, cell, subsys=subsys)
1822 CALL xyz_outerdiag_colvar(colvar, cell, subsys=subsys)
1824 CALL u_colvar(colvar, force_env=force_env)
1826 CALL wc_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1828 CALL hbp_colvar(colvar, cell, subsys=subsys, qs_env=qs_env)
1830 CALL ring_puckering_colvar(colvar, cell, subsys=subsys)
1832 CALL mindist_colvar(colvar, cell, subsys=subsys)
1834 CALL colvar_eval_voronoiipz(colvar, cell, subsys=subsys)
1836 cpabort(
"Unknown colvar type for colvar_eval_glob_f")
1850 SUBROUTINE colvar_recursive_eval(colvar, cell, particles)
1857 colvar%dsdr = 0.0_dp
1858 SELECT CASE (colvar%type_id)
1860 CALL dist_colvar(colvar, cell, particles=particles)
1862 CALL coord_colvar(colvar, cell, particles=particles)
1864 CALL torsion_colvar(colvar, cell, particles=particles)
1866 CALL angle_colvar(colvar, cell, particles=particles)
1868 CALL dfunct_colvar(colvar, cell, particles=particles)
1870 CALL plane_distance_colvar(colvar, cell, particles=particles)
1872 CALL plane_plane_angle_colvar(colvar, cell, particles=particles)
1874 CALL rotation_colvar(colvar, cell, particles=particles)
1876 CALL qparm_colvar(colvar, cell, particles=particles)
1878 CALL hydronium_shell_colvar(colvar, cell, particles=particles)
1880 CALL hydronium_dist_colvar(colvar, cell, particles=particles)
1882 CALL acid_hyd_dist_colvar(colvar, cell, particles=particles)
1884 CALL acid_hyd_shell_colvar(colvar, cell, particles=particles)
1886 CALL rmsd_colvar(colvar, particles=particles)
1888 CALL reaction_path_colvar(colvar, cell, particles=particles)
1890 CALL distance_from_path_colvar(colvar, cell, particles=particles)
1892 CALL combine_colvar(colvar, cell, particles=particles)
1894 CALL xyz_diag_colvar(colvar, cell, particles=particles)
1896 CALL xyz_outerdiag_colvar(colvar, cell, particles=particles)
1898 CALL ring_puckering_colvar(colvar, cell, particles=particles)
1900 CALL mindist_colvar(colvar, cell, particles=particles)
1902 cpabort(
"need force_env!")
1904 CALL wc_colvar(colvar, cell, particles=particles)
1906 CALL hbp_colvar(colvar, cell, particles=particles)
1908 CALL colvar_eval_voronoiipz(colvar, cell, particles=particles)
1910 cpabort(
"Unknown colvar type for colvar_recursive_eval")
1912 END SUBROUTINE colvar_recursive_eval
1922 SUBROUTINE get_coordinates(colvar, i, ri, my_particles)
1924 INTEGER,
INTENT(IN) :: i
1925 REAL(kind=
dp),
DIMENSION(3),
INTENT(OUT) :: ri
1928 IF (colvar%use_points)
THEN
1931 ri(:) = my_particles(i)%r(:)
1934 END SUBROUTINE get_coordinates
1944 SUBROUTINE get_mass(colvar, i, mi, my_particles)
1946 INTEGER,
INTENT(IN) :: i
1947 REAL(kind=
dp),
INTENT(OUT) :: mi
1950 IF (colvar%use_points)
THEN
1953 mi = my_particles(i)%atomic_kind%mass
1956 END SUBROUTINE get_mass
1965 SUBROUTINE put_derivative(colvar, i, fi)
1967 INTEGER,
INTENT(IN) :: i
1968 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: fi
1970 IF (colvar%use_points)
THEN
1973 colvar%dsdr(:, i) = colvar%dsdr(:, i) + fi
1976 END SUBROUTINE put_derivative
1986 SUBROUTINE xyz_diag_colvar(colvar, cell, subsys, particles)
1991 POINTER :: particles
1994 REAL(
dp) :: fi(3), r, r0(3), ss(3), xi(3), xpi(3)
1998 NULLIFY (particles_i)
2001 IF (
PRESENT(particles))
THEN
2002 my_particles => particles
2004 cpassert(
PRESENT(subsys))
2006 my_particles => particles_i%els
2008 i = colvar%xyz_diag_param%i_atom
2010 CALL get_coordinates(colvar, i, xpi, my_particles)
2013 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
2014 IF (all(colvar%xyz_diag_param%r0 == huge(0.0_dp)))
THEN
2015 colvar%xyz_diag_param%r0 = xpi
2017 r0 = colvar%xyz_diag_param%r0
2022 IF (colvar%xyz_diag_param%use_pbc)
THEN
2023 ss = matmul(cell%h_inv, xpi - r0)
2025 xi = matmul(cell%hmat, ss)
2030 IF (.NOT. colvar%xyz_diag_param%use_absolute_position)
THEN
2031 SELECT CASE (colvar%xyz_diag_param%component)
2051 r = xi(1)**2 + xi(2)**2 + xi(3)**2
2054 SELECT CASE (colvar%xyz_diag_param%component)
2071 cpabort(
"xyz_diag_colvar not implemented for anything which is not a single component")
2077 CALL put_derivative(colvar, 1, fi)
2079 END SUBROUTINE xyz_diag_colvar
2089 SUBROUTINE xyz_outerdiag_colvar(colvar, cell, subsys, particles)
2094 POINTER :: particles
2097 REAL(
dp) :: fi(3, 2), r, r0(3), ss(3), xi(3, 2), &
2102 NULLIFY (particles_i)
2105 IF (
PRESENT(particles))
THEN
2106 my_particles => particles
2108 cpassert(
PRESENT(subsys))
2110 my_particles => particles_i%els
2113 i = colvar%xyz_outerdiag_param%i_atoms(k)
2115 CALL get_coordinates(colvar, i, xpi, my_particles)
2116 r0 = colvar%xyz_outerdiag_param%r0(:, k)
2117 IF (all(colvar%xyz_outerdiag_param%r0(:, k) == huge(0.0_dp))) r0 = xpi
2119 IF (colvar%xyz_outerdiag_param%use_pbc)
THEN
2120 ss = matmul(cell%h_inv, xpi - r0)
2122 xi(:, k) = matmul(cell%hmat, ss)
2127 SELECT CASE (colvar%xyz_outerdiag_param%components(k))
2152 IF (xi(l, 1) /= 0.0_dp) fi(l, 1) = fi(l, 1) + xi(i, 2)
2153 r = r + xi(l, 1)*xi(i, 2)
2155 IF (xi(i, 2) /= 0.0_dp) fi(i, 2) = sum(xi(:, 1))
2159 CALL put_derivative(colvar, 1, fi(:, 1))
2160 CALL put_derivative(colvar, 2, fi(:, 2))
2162 END SUBROUTINE xyz_outerdiag_colvar
2172 SUBROUTINE u_colvar(colvar, force_env)
2176 CHARACTER(LEN=default_path_length) :: coupling_function
2177 CHARACTER(LEN=default_string_length) :: def_error, this_error
2178 CHARACTER(LEN=default_string_length), &
2179 ALLOCATABLE,
DIMENSION(:) :: parameters
2180 INTEGER :: iatom, iforce_eval, iparticle, &
2181 jparticle, natom, natom_iforce, &
2183 INTEGER,
DIMENSION(:),
POINTER :: glob_natoms, map_index
2184 REAL(
dp) :: dedf, dx, err, fi(3), lerr, &
2186 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: values
2195 NULLIFY (particles_main, subsys_main)
2197 CALL cp_subsys_get(subsys=subsys_main, particles=particles_main)
2198 natom =
SIZE(particles_main%els)
2199 colvar%n_atom_s = natom
2200 colvar%u_param%natom = natom
2204 colvar%i_atom(iatom) = iatom
2207 IF (.NOT.
ASSOCIATED(colvar%u_param%mixed_energy_section))
THEN
2208 CALL force_env_get(force_env, potential_energy=potential_energy)
2209 colvar%ss = potential_energy
2213 fi(:) = -particles_main%els(iatom)%f
2214 CALL put_derivative(colvar, iatom, fi)
2218 CALL cp_abort(__location__, &
2219 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
2220 ' A combination of mixed force_eval energies has been requested as '// &
2221 ' collective variable, but the MIXED env is not in use! Aborting.')
2223 CALL force_env_get(force_env, force_env_section=force_env_section)
2225 NULLIFY (subsystems, particles, global_forces, map_index, glob_natoms)
2226 nforce_eval =
SIZE(force_env%sub_force_env)
2227 ALLOCATE (glob_natoms(nforce_eval))
2228 ALLOCATE (subsystems(nforce_eval))
2229 ALLOCATE (particles(nforce_eval))
2231 ALLOCATE (global_forces(nforce_eval))
2234 DO iforce_eval = 1, nforce_eval
2235 NULLIFY (subsystems(iforce_eval)%subsys, particles(iforce_eval)%list)
2236 IF (.NOT.
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
2238 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env, &
2239 subsys=subsystems(iforce_eval)%subsys)
2242 particles=particles(iforce_eval)%list)
2245 natom_iforce =
SIZE(particles(iforce_eval)%list%els)
2248 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2249 glob_natoms(iforce_eval) = natom_iforce
2254 CALL force_env%para_env%sync()
2255 CALL force_env%para_env%sum(glob_natoms)
2258 DO iforce_eval = 1, nforce_eval
2259 ALLOCATE (global_forces(iforce_eval)%forces(3, glob_natoms(iforce_eval)))
2260 global_forces(iforce_eval)%forces = 0.0_dp
2261 IF (
ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env))
THEN
2262 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source())
THEN
2264 DO iparticle = 1, glob_natoms(iforce_eval)
2265 global_forces(iforce_eval)%forces(:, iparticle) = &
2266 particles(iforce_eval)%list%els(iparticle)%f
2270 CALL force_env%para_env%sum(global_forces(iforce_eval)%forces)
2273 wrk_section => colvar%u_param%mixed_energy_section
2276 parameters, values, var_values=force_env%mixed_env%energies, &
2277 dx=dx, error_limit=lerr)
2279 CALL parsef(1, trim(coupling_function), parameters)
2281 colvar%ss =
evalf(1, values)
2284 DO iforce_eval = 1, nforce_eval
2285 dedf =
evalfd(1, iforce_eval, values, dx, err)
2286 IF (abs(err) > lerr)
THEN
2287 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
2288 WRITE (def_error,
"(A,G12.6,A)")
"(", lerr,
")"
2291 CALL cp_warn(__location__, &
2292 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
2293 ' Error '//trim(this_error)//
' in computing numerical derivatives larger than'// &
2294 trim(def_error)//
' .')
2299 nforce_eval, map_index)
2302 DO iparticle = 1, glob_natoms(iforce_eval)
2303 jparticle = map_index(iparticle)
2304 fi = -dedf*global_forces(iforce_eval)%forces(:, iparticle)
2305 CALL put_derivative(colvar, jparticle, fi)
2308 IF (
ASSOCIATED(map_index))
THEN
2309 DEALLOCATE (map_index)
2313 DO iforce_eval = 1, nforce_eval
2314 DEALLOCATE (global_forces(iforce_eval)%forces)
2316 DEALLOCATE (glob_natoms)
2318 DEALLOCATE (parameters)
2319 DEALLOCATE (global_forces)
2320 DEALLOCATE (subsystems)
2321 DEALLOCATE (particles)
2324 END SUBROUTINE u_colvar
2334 SUBROUTINE plane_distance_colvar(colvar, cell, subsys, particles)
2340 POINTER :: particles
2342 INTEGER :: i, j, k, l
2343 REAL(
dp) :: a, b, dsdxpn(3), dxpndxi(3, 3), dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), &
2344 fk(3), fl(3), r12, ri(3), rj(3), rk(3), rl(3), ss(3), xpij(3), xpkj(3), xpl(3), xpn(3)
2348 NULLIFY (particles_i)
2351 IF (
PRESENT(particles))
THEN
2352 my_particles => particles
2354 cpassert(
PRESENT(subsys))
2356 my_particles => particles_i%els
2358 i = colvar%plane_distance_param%plane(1)
2359 j = colvar%plane_distance_param%plane(2)
2360 k = colvar%plane_distance_param%plane(3)
2361 l = colvar%plane_distance_param%point
2363 CALL get_coordinates(colvar, i, ri, my_particles)
2364 CALL get_coordinates(colvar, j, rj, my_particles)
2365 CALL get_coordinates(colvar, k, rk, my_particles)
2366 CALL get_coordinates(colvar, l, rl, my_particles)
2369 xpl = rl - (ri + rj + rk)/3.0_dp
2370 IF (colvar%plane_distance_param%use_pbc)
THEN
2372 ss = matmul(cell%h_inv, ri - rj)
2374 xpij = matmul(cell%hmat, ss)
2376 ss = matmul(cell%h_inv, rk - rj)
2378 xpkj = matmul(cell%hmat, ss)
2380 ss = matmul(cell%h_inv, rl - (ri + rj + rk)/3.0_dp)
2382 xpl = matmul(cell%hmat, ss)
2385 xpn(1) = xpij(2)*xpkj(3) - xpij(3)*xpkj(2)
2386 xpn(2) = xpij(3)*xpkj(1) - xpij(1)*xpkj(3)
2387 xpn(3) = xpij(1)*xpkj(2) - xpij(2)*xpkj(1)
2388 a = dot_product(xpn, xpn)
2389 b = dot_product(xpl, xpn)
2392 dsdxpn(1) = xpl(1)/r12 - b*xpn(1)/(r12*a)
2393 dsdxpn(2) = xpl(2)/r12 - b*xpn(2)/(r12*a)
2394 dsdxpn(3) = xpl(3)/r12 - b*xpn(3)/(r12*a)
2396 dxpndxi(1, 1) = 0.0_dp
2397 dxpndxi(1, 2) = 1.0_dp*xpkj(3)
2398 dxpndxi(1, 3) = -1.0_dp*xpkj(2)
2399 dxpndxi(2, 1) = -1.0_dp*xpkj(3)
2400 dxpndxi(2, 2) = 0.0_dp
2401 dxpndxi(2, 3) = 1.0_dp*xpkj(1)
2402 dxpndxi(3, 1) = 1.0_dp*xpkj(2)
2403 dxpndxi(3, 2) = -1.0_dp*xpkj(1)
2404 dxpndxi(3, 3) = 0.0_dp
2406 dxpndxj(1, 1) = 0.0_dp
2407 dxpndxj(1, 2) = -1.0_dp*xpkj(3) + xpij(3)
2408 dxpndxj(1, 3) = -1.0_dp*xpij(2) + xpkj(2)
2409 dxpndxj(2, 1) = -1.0_dp*xpij(3) + xpkj(3)
2410 dxpndxj(2, 2) = 0.0_dp
2411 dxpndxj(2, 3) = -1.0_dp*xpkj(1) + xpij(1)
2412 dxpndxj(3, 1) = -1.0_dp*xpkj(2) + xpij(2)
2413 dxpndxj(3, 2) = -1.0_dp*xpij(1) + xpkj(1)
2414 dxpndxj(3, 3) = 0.0_dp
2416 dxpndxk(1, 1) = 0.0_dp
2417 dxpndxk(1, 2) = -1.0_dp*xpij(3)
2418 dxpndxk(1, 3) = 1.0_dp*xpij(2)
2419 dxpndxk(2, 1) = 1.0_dp*xpij(3)
2420 dxpndxk(2, 2) = 0.0_dp
2421 dxpndxk(2, 3) = -1.0_dp*xpij(1)
2422 dxpndxk(3, 1) = -1.0_dp*xpij(2)
2423 dxpndxk(3, 2) = 1.0_dp*xpij(1)
2424 dxpndxk(3, 3) = 0.0_dp
2426 fi(:) = matmul(dsdxpn, dxpndxi) - xpn/(3.0_dp*r12)
2427 fj(:) = matmul(dsdxpn, dxpndxj) - xpn/(3.0_dp*r12)
2428 fk(:) = matmul(dsdxpn, dxpndxk) - xpn/(3.0_dp*r12)
2431 CALL put_derivative(colvar, 1, fi)
2432 CALL put_derivative(colvar, 2, fj)
2433 CALL put_derivative(colvar, 3, fk)
2434 CALL put_derivative(colvar, 4, fl)
2436 END SUBROUTINE plane_distance_colvar
2447 SUBROUTINE plane_plane_angle_colvar(colvar, cell, subsys, particles)
2453 POINTER :: particles
2455 INTEGER :: i1, i2, j1, j2, k1, k2, np
2457 REAL(
dp) :: a1, a2, d, dnorm_dxpn(3), dprod12_dxpn(3), dsdxpn(3), dt_dxpn(3), dxpndxi(3, 3), &
2458 dxpndxj(3, 3), dxpndxk(3, 3), fi(3), fj(3), fk(3), fmod, norm1, norm2, prod_12, ri1(3), &
2459 ri2(3), rj1(3), rj2(3), rk1(3), rk2(3), ss(3), t, xpij1(3), xpij2(3), xpkj1(3), xpkj2(3), &
2464 NULLIFY (particles_i)
2468 IF (
PRESENT(particles))
THEN
2469 my_particles => particles
2471 cpassert(
PRESENT(subsys))
2473 my_particles => particles_i%els
2477 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2478 i1 = colvar%plane_plane_angle_param%plane1%points(1)
2479 j1 = colvar%plane_plane_angle_param%plane1%points(2)
2480 k1 = colvar%plane_plane_angle_param%plane1%points(3)
2483 CALL get_coordinates(colvar, i1, ri1, my_particles)
2484 CALL get_coordinates(colvar, j1, rj1, my_particles)
2485 CALL get_coordinates(colvar, k1, rk1, my_particles)
2488 ss = matmul(cell%h_inv, ri1 - rj1)
2490 xpij1 = matmul(cell%hmat, ss)
2493 ss = matmul(cell%h_inv, rk1 - rj1)
2495 xpkj1 = matmul(cell%hmat, ss)
2498 xpn1(1) = xpij1(2)*xpkj1(3) - xpij1(3)*xpkj1(2)
2499 xpn1(2) = xpij1(3)*xpkj1(1) - xpij1(1)*xpkj1(3)
2500 xpn1(3) = xpij1(1)*xpkj1(2) - xpij1(2)*xpkj1(1)
2502 xpn1 = colvar%plane_plane_angle_param%plane1%normal_vec
2504 a1 = dot_product(xpn1, xpn1)
2506 cpassert(norm1 /= 0.0_dp)
2509 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2510 i2 = colvar%plane_plane_angle_param%plane2%points(1)
2511 j2 = colvar%plane_plane_angle_param%plane2%points(2)
2512 k2 = colvar%plane_plane_angle_param%plane2%points(3)
2515 CALL get_coordinates(colvar, i2, ri2, my_particles)
2516 CALL get_coordinates(colvar, j2, rj2, my_particles)
2517 CALL get_coordinates(colvar, k2, rk2, my_particles)
2520 ss = matmul(cell%h_inv, ri2 - rj2)
2522 xpij2 = matmul(cell%hmat, ss)
2525 ss = matmul(cell%h_inv, rk2 - rj2)
2527 xpkj2 = matmul(cell%hmat, ss)
2530 xpn2(1) = xpij2(2)*xpkj2(3) - xpij2(3)*xpkj2(2)
2531 xpn2(2) = xpij2(3)*xpkj2(1) - xpij2(1)*xpkj2(3)
2532 xpn2(3) = xpij2(1)*xpkj2(2) - xpij2(2)*xpkj2(1)
2534 xpn2 = colvar%plane_plane_angle_param%plane2%normal_vec
2536 a2 = dot_product(xpn2, xpn2)
2538 cpassert(norm2 /= 0.0_dp)
2541 prod_12 = dot_product(xpn1, xpn2)
2545 t = min(1.0_dp, abs(t))*sign(1.0_dp, t)
2548 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2551 fmod = -1.0_dp/sin(colvar%ss)
2556 IF (colvar%plane_plane_angle_param%plane1%type_of_def ==
plane_def_atoms)
THEN
2558 dnorm_dxpn = 1.0_dp/norm1*xpn1
2559 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm2)/d**2
2561 dsdxpn(1) = fmod*dt_dxpn(1)
2562 dsdxpn(2) = fmod*dt_dxpn(2)
2563 dsdxpn(3) = fmod*dt_dxpn(3)
2565 dxpndxi(1, 1) = 0.0_dp
2566 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2567 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2568 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2569 dxpndxi(2, 2) = 0.0_dp
2570 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2571 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2572 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2573 dxpndxi(3, 3) = 0.0_dp
2575 dxpndxj(1, 1) = 0.0_dp
2576 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2577 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2578 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2579 dxpndxj(2, 2) = 0.0_dp
2580 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2581 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2582 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2583 dxpndxj(3, 3) = 0.0_dp
2585 dxpndxk(1, 1) = 0.0_dp
2586 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2587 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2588 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2589 dxpndxk(2, 2) = 0.0_dp
2590 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2591 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2592 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2593 dxpndxk(3, 3) = 0.0_dp
2595 fi = matmul(dsdxpn, dxpndxi)
2596 fj = matmul(dsdxpn, dxpndxj)
2597 fk = matmul(dsdxpn, dxpndxk)
2600 CALL put_derivative(colvar, np + 1, fi)
2601 CALL put_derivative(colvar, np + 2, fj)
2602 CALL put_derivative(colvar, np + 3, fk)
2607 IF (colvar%plane_plane_angle_param%plane2%type_of_def ==
plane_def_atoms)
THEN
2609 dnorm_dxpn = 1.0_dp/norm2*xpn2
2610 dt_dxpn = (dprod12_dxpn*d - prod_12*dnorm_dxpn*norm1)/d**2
2612 dsdxpn(1) = fmod*dt_dxpn(1)
2613 dsdxpn(2) = fmod*dt_dxpn(2)
2614 dsdxpn(3) = fmod*dt_dxpn(3)
2616 dxpndxi(1, 1) = 0.0_dp
2617 dxpndxi(1, 2) = 1.0_dp*xpkj1(3)
2618 dxpndxi(1, 3) = -1.0_dp*xpkj1(2)
2619 dxpndxi(2, 1) = -1.0_dp*xpkj1(3)
2620 dxpndxi(2, 2) = 0.0_dp
2621 dxpndxi(2, 3) = 1.0_dp*xpkj1(1)
2622 dxpndxi(3, 1) = 1.0_dp*xpkj1(2)
2623 dxpndxi(3, 2) = -1.0_dp*xpkj1(1)
2624 dxpndxi(3, 3) = 0.0_dp
2626 dxpndxj(1, 1) = 0.0_dp
2627 dxpndxj(1, 2) = -1.0_dp*xpkj1(3) + xpij1(3)
2628 dxpndxj(1, 3) = -1.0_dp*xpij1(2) + xpkj1(2)
2629 dxpndxj(2, 1) = -1.0_dp*xpij1(3) + xpkj1(3)
2630 dxpndxj(2, 2) = 0.0_dp
2631 dxpndxj(2, 3) = -1.0_dp*xpkj1(1) + xpij1(1)
2632 dxpndxj(3, 1) = -1.0_dp*xpkj1(2) + xpij1(2)
2633 dxpndxj(3, 2) = -1.0_dp*xpij1(1) + xpkj1(1)
2634 dxpndxj(3, 3) = 0.0_dp
2636 dxpndxk(1, 1) = 0.0_dp
2637 dxpndxk(1, 2) = -1.0_dp*xpij1(3)
2638 dxpndxk(1, 3) = 1.0_dp*xpij1(2)
2639 dxpndxk(2, 1) = 1.0_dp*xpij1(3)
2640 dxpndxk(2, 2) = 0.0_dp
2641 dxpndxk(2, 3) = -1.0_dp*xpij1(1)
2642 dxpndxk(3, 1) = -1.0_dp*xpij1(2)
2643 dxpndxk(3, 2) = 1.0_dp*xpij1(1)
2644 dxpndxk(3, 3) = 0.0_dp
2646 fi = matmul(dsdxpn, dxpndxi)
2647 fj = matmul(dsdxpn, dxpndxj)
2648 fk = matmul(dsdxpn, dxpndxk)
2651 CALL put_derivative(colvar, np + 1, fi)
2652 CALL put_derivative(colvar, np + 2, fj)
2653 CALL put_derivative(colvar, np + 3, fk)
2656 END SUBROUTINE plane_plane_angle_colvar
2666 SUBROUTINE rotation_colvar(colvar, cell, subsys, particles)
2671 POINTER :: particles
2674 REAL(
dp) :: a, b, fmod, t0, t1, t2, t3, xdum(3), &
2676 REAL(kind=
dp) :: dp1b1(3), dp1b2(3), dp2b1(3), dp2b2(3), &
2677 ss(3), xp1b1(3), xp1b2(3), xp2b1(3), &
2682 NULLIFY (particles_i)
2685 IF (
PRESENT(particles))
THEN
2686 my_particles => particles
2688 cpassert(
PRESENT(subsys))
2690 my_particles => particles_i%els
2692 i = colvar%rotation_param%i_at1_bond1
2693 CALL get_coordinates(colvar, i, xp1b1, my_particles)
2694 i = colvar%rotation_param%i_at2_bond1
2695 CALL get_coordinates(colvar, i, xp2b1, my_particles)
2696 i = colvar%rotation_param%i_at1_bond2
2697 CALL get_coordinates(colvar, i, xp1b2, my_particles)
2698 i = colvar%rotation_param%i_at2_bond2
2699 CALL get_coordinates(colvar, i, xp2b2, my_particles)
2701 ss = matmul(cell%h_inv, xp1b1 - xp2b1)
2703 xij = matmul(cell%hmat, ss)
2705 ss = matmul(cell%h_inv, xp1b2 - xp2b2)
2707 xkj = matmul(cell%hmat, ss)
2712 t1 = 1.0_dp/(a**3.0_dp*b)
2713 t2 = 1.0_dp/(a*b**3.0_dp)
2714 t3 = dot_product(xij, xkj)
2715 colvar%ss = acos(t3*t0)
2716 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2719 fmod = -1.0_dp/sin(colvar%ss)
2721 dp1b1 = xkj(:)*t0 - xij(:)*t1*t3
2722 dp2b1 = -xkj(:)*t0 + xij(:)*t1*t3
2723 dp1b2 = xij(:)*t0 - xkj(:)*t2*t3
2724 dp2b2 = -xij(:)*t0 + xkj(:)*t2*t3
2727 idum = colvar%rotation_param%i_at1_bond1
2728 CALL put_derivative(colvar, idum, xdum)
2730 idum = colvar%rotation_param%i_at2_bond1
2731 CALL put_derivative(colvar, idum, xdum)
2733 idum = colvar%rotation_param%i_at1_bond2
2734 CALL put_derivative(colvar, idum, xdum)
2736 idum = colvar%rotation_param%i_at2_bond2
2737 CALL put_derivative(colvar, idum, xdum)
2739 END SUBROUTINE rotation_colvar
2750 SUBROUTINE dfunct_colvar(colvar, cell, subsys, particles)
2755 POINTER :: particles
2757 INTEGER :: i, j, k, l
2758 REAL(
dp) :: fi(3), fj(3), fk(3), fl(3), r12, r34, &
2759 ss(3), xij(3), xkl(3), xpi(3), xpj(3), &
2764 NULLIFY (particles_i)
2767 IF (
PRESENT(particles))
THEN
2768 my_particles => particles
2770 cpassert(
PRESENT(subsys))
2772 my_particles => particles_i%els
2774 i = colvar%dfunct_param%i_at_dfunct(1)
2775 j = colvar%dfunct_param%i_at_dfunct(2)
2777 CALL get_coordinates(colvar, i, xpi, my_particles)
2778 CALL get_coordinates(colvar, j, xpj, my_particles)
2779 IF (colvar%dfunct_param%use_pbc)
THEN
2780 ss = matmul(cell%h_inv, xpi - xpj)
2782 xij = matmul(cell%hmat, ss)
2786 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2788 k = colvar%dfunct_param%i_at_dfunct(3)
2789 l = colvar%dfunct_param%i_at_dfunct(4)
2790 CALL get_coordinates(colvar, k, xpk, my_particles)
2791 CALL get_coordinates(colvar, l, xpl, my_particles)
2792 IF (colvar%dfunct_param%use_pbc)
THEN
2793 ss = matmul(cell%h_inv, xpk - xpl)
2795 xkl = matmul(cell%hmat, ss)
2799 r34 = sqrt(xkl(1)**2 + xkl(2)**2 + xkl(3)**2)
2801 colvar%ss = r12 + colvar%dfunct_param%coeff*r34
2804 fk(:) = colvar%dfunct_param%coeff*xkl/r34
2805 fl(:) = -colvar%dfunct_param%coeff*xkl/r34
2806 CALL put_derivative(colvar, 1, fi)
2807 CALL put_derivative(colvar, 2, fj)
2808 CALL put_derivative(colvar, 3, fk)
2809 CALL put_derivative(colvar, 4, fl)
2811 END SUBROUTINE dfunct_colvar
2821 SUBROUTINE angle_colvar(colvar, cell, subsys, particles)
2826 POINTER :: particles
2829 REAL(
dp) :: a, b, fi(3), fj(3), fk(3), fmod, ri(3), &
2830 rj(3), rk(3), ss(3), t0, t1, t2, t3, &
2835 NULLIFY (particles_i)
2838 IF (
PRESENT(particles))
THEN
2839 my_particles => particles
2841 cpassert(
PRESENT(subsys))
2843 my_particles => particles_i%els
2845 i = colvar%angle_param%i_at_angle(1)
2846 j = colvar%angle_param%i_at_angle(2)
2847 k = colvar%angle_param%i_at_angle(3)
2848 CALL get_coordinates(colvar, i, ri, my_particles)
2849 CALL get_coordinates(colvar, j, rj, my_particles)
2850 CALL get_coordinates(colvar, k, rk, my_particles)
2852 ss = matmul(cell%h_inv, ri - rj)
2854 xij = matmul(cell%hmat, ss)
2856 ss = matmul(cell%h_inv, rk - rj)
2858 xkj = matmul(cell%hmat, ss)
2863 t1 = 1.0_dp/(a**3.0_dp*b)
2864 t2 = 1.0_dp/(a*b**3.0_dp)
2865 t3 = dot_product(xij, xkj)
2866 colvar%ss = acos(t3*t0)
2867 IF ((abs(colvar%ss) < tolerance_acos) .OR. (abs(colvar%ss -
pi) < tolerance_acos))
THEN
2870 fmod = -1.0_dp/sin(colvar%ss)
2872 fi(:) = xkj(:)*t0 - xij(:)*t1*t3
2873 fj(:) = -xkj(:)*t0 + xij(:)*t1*t3 - xij(:)*t0 + xkj(:)*t2*t3
2874 fk(:) = xij(:)*t0 - xkj(:)*t2*t3
2878 CALL put_derivative(colvar, 1, fi)
2879 CALL put_derivative(colvar, 2, fj)
2880 CALL put_derivative(colvar, 3, fk)
2882 END SUBROUTINE angle_colvar
2892 SUBROUTINE dist_colvar(colvar, cell, subsys, particles)
2897 POINTER :: particles
2900 REAL(
dp) :: fi(3), fj(3), r12, ss(3), xij(3), &
2905 NULLIFY (particles_i)
2908 IF (
PRESENT(particles))
THEN
2909 my_particles => particles
2911 cpassert(
PRESENT(subsys))
2913 my_particles => particles_i%els
2915 i = colvar%dist_param%i_at
2916 j = colvar%dist_param%j_at
2917 CALL get_coordinates(colvar, i, xpi, my_particles)
2918 CALL get_coordinates(colvar, j, xpj, my_particles)
2919 ss = matmul(cell%h_inv, xpi - xpj)
2921 xij = matmul(cell%hmat, ss)
2922 SELECT CASE (colvar%dist_param%axis_id)
2941 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
2943 IF (colvar%dist_param%sign_d)
THEN
2944 SELECT CASE (colvar%dist_param%axis_id)
2962 CALL put_derivative(colvar, 1, fi)
2963 CALL put_derivative(colvar, 2, fj)
2965 END SUBROUTINE dist_colvar
2976 SUBROUTINE torsion_colvar(colvar, cell, subsys, particles, no_riemann_sheet_op)
2982 POINTER :: particles
2983 LOGICAL,
INTENT(IN),
OPTIONAL :: no_riemann_sheet_op
2986 LOGICAL :: no_riemann_sheet
2987 REAL(
dp) ::
angle, cosine, dedphi, dedxia, dedxib, dedxic, dedxid, dedxt, dedxu, dedyia, &
2988 dedyib, dedyic, dedyid, dedyt, dedyu, dedzia, dedzib, dedzic, dedzid, dedzt, dedzu, dt, &
2989 e, ftmp(3), o0, rcb, rt2, rtmp(3), rtru, ru2, sine, ss(3), xba, xca, xcb, xdb, xdc, xt, &
2990 xtu, xu, yba, yca, ycb, ydb, ydc, yt, ytu, yu, zba, zca, zcb, zdb, zdc, zt, ztu, zu
2991 REAL(
dp),
DIMENSION(3, 4) :: rr
2995 NULLIFY (particles_i)
2997 IF (
PRESENT(particles))
THEN
2998 my_particles => particles
3000 cpassert(
PRESENT(subsys))
3002 my_particles => particles_i%els
3004 no_riemann_sheet = .false.
3005 IF (
PRESENT(no_riemann_sheet_op)) no_riemann_sheet = no_riemann_sheet_op
3007 i = colvar%torsion_param%i_at_tors(ii)
3008 CALL get_coordinates(colvar, i, rtmp, my_particles)
3009 rr(:, ii) = rtmp(1:3)
3011 o0 = colvar%torsion_param%o0
3013 ss = matmul(cell%h_inv, rr(:, 2) - rr(:, 1))
3015 ss = matmul(cell%hmat, ss)
3020 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 2))
3022 ss = matmul(cell%hmat, ss)
3027 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 3))
3029 ss = matmul(cell%hmat, ss)
3034 xt = yba*zcb - ycb*zba
3035 yt = zba*xcb - zcb*xba
3036 zt = xba*ycb - xcb*yba
3037 xu = ycb*zdc - ydc*zcb
3038 yu = zcb*xdc - zdc*xcb
3039 zu = xcb*ydc - xdc*ycb
3043 rt2 = xt*xt + yt*yt + zt*zt
3044 ru2 = xu*xu + yu*yu + zu*zu
3045 rtru = sqrt(rt2*ru2)
3046 IF (rtru /= 0.0_dp)
THEN
3047 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb)
3048 cosine = (xt*xu + yt*yu + zt*zu)/rtru
3049 sine = (xcb*xtu + ycb*ytu + zcb*ztu)/(rcb*rtru)
3050 cosine = min(1.0_dp, max(-1.0_dp, cosine))
3051 angle = acos(cosine)
3055 dt = mod(2.0e4_dp*
pi + dt - o0, 2.0_dp*
pi)
3056 IF (dt >
pi) dt = dt - 2.0_dp*
pi
3058 colvar%torsion_param%o0 = dt
3068 ss = matmul(cell%h_inv, rr(:, 3) - rr(:, 1))
3070 ss = matmul(cell%hmat, ss)
3075 ss = matmul(cell%h_inv, rr(:, 4) - rr(:, 2))
3077 ss = matmul(cell%hmat, ss)
3082 dedxt = dedphi*(yt*zcb - ycb*zt)/(rt2*rcb)
3083 dedyt = dedphi*(zt*xcb - zcb*xt)/(rt2*rcb)
3084 dedzt = dedphi*(xt*ycb - xcb*yt)/(rt2*rcb)
3085 dedxu = -dedphi*(yu*zcb - ycb*zu)/(ru2*rcb)
3086 dedyu = -dedphi*(zu*xcb - zcb*xu)/(ru2*rcb)
3087 dedzu = -dedphi*(xu*ycb - xcb*yu)/(ru2*rcb)
3091 dedxia = zcb*dedyt - ycb*dedzt
3092 dedyia = xcb*dedzt - zcb*dedxt
3093 dedzia = ycb*dedxt - xcb*dedyt
3094 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu
3095 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu
3096 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu
3097 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu
3098 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu
3099 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu
3100 dedxid = zcb*dedyu - ycb*dedzu
3101 dedyid = xcb*dedzu - zcb*dedxu
3102 dedzid = ycb*dedxu - xcb*dedyu
3119 IF (no_riemann_sheet) colvar%ss = atan2(sin(e), cos(e))
3123 CALL put_derivative(colvar, 1, ftmp)
3127 CALL put_derivative(colvar, 2, ftmp)
3131 CALL put_derivative(colvar, 3, ftmp)
3135 CALL put_derivative(colvar, 4, ftmp)
3136 END SUBROUTINE torsion_colvar
3145 SUBROUTINE qparm_colvar(colvar, cell, subsys, particles)
3150 POINTER :: particles
3152 INTEGER :: aa, bb, cc, i, idim, ii, j, jj, l, mm, &
3153 n_atoms_from, n_atoms_to, ncells(3)
3154 LOGICAL :: include_images
3155 REAL(kind=
dp) :: denominator_tolerance, fact, ftmp(3), im_qlm, inv_n_atoms_from, nbond, &
3156 pre_fac, ql, qparm, r1cut, rcut, re_qlm, rij, rij_shift, shift(3), ss(3), ss0(3), xij(3), &
3158 REAL(kind=
dp),
DIMENSION(3) :: d_im_qlm_dxi, d_nbond_dxi, d_ql_dxi, &
3159 d_re_qlm_dxi, xpi, xpj
3163 n_atoms_to = colvar%qparm_param%n_atoms_to
3164 n_atoms_from = colvar%qparm_param%n_atoms_from
3165 rcut = colvar%qparm_param%rcut
3166 l = colvar%qparm_param%l
3167 r1cut = colvar%qparm_param%rstart
3168 include_images = colvar%qparm_param%include_images
3169 NULLIFY (particles_i)
3171 IF (
PRESENT(particles))
THEN
3172 my_particles => particles
3174 cpassert(
PRESENT(subsys))
3176 my_particles => particles_i%els
3178 cpassert(r1cut < rcut)
3179 denominator_tolerance = 1.0e-8_dp
3182 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
3183 DO ii = 1, n_atoms_from
3184 i = colvar%qparm_param%i_at_from(ii)
3185 CALL get_coordinates(colvar, i, xpi, my_particles)
3188 d_ql_dxi(:) = 0.0_dp
3194 d_re_qlm_dxi(:) = 0.0_dp
3195 d_im_qlm_dxi(:) = 0.0_dp
3196 d_nbond_dxi(:) = 0.0_dp
3198 jloop:
DO jj = 1, n_atoms_to
3200 j = colvar%qparm_param%i_at_to(jj)
3201 CALL get_coordinates(colvar, j, xpj, my_particles)
3203 IF (include_images)
THEN
3205 cpassert(cell%orthorhombic)
3209 xij(:) = xpj(:) - xpi(:)
3210 ss = matmul(cell%h_inv, xij)
3216 shift(idim) = 1.0_dp
3217 xij_shift = matmul(cell%hmat, shift)
3218 rij_shift = norm2(xij_shift)
3219 ncells(idim) = floor(rcut/rij_shift - 0.5_dp)
3223 DO aa = -ncells(1), ncells(1)
3224 DO bb = -ncells(2), ncells(2)
3225 DO cc = -ncells(3), ncells(3)
3227 IF (i == j .AND. aa == 0 .AND. bb == 0 .AND. cc == 0) cycle
3228 shift(1) = real(aa, kind=
dp)
3229 shift(2) = real(bb, kind=
dp)
3230 shift(3) = real(cc, kind=
dp)
3231 xij = matmul(cell%hmat, ss0(:) + shift(:))
3233 IF (rij > rcut) cycle
3236 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3237 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3238 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3246 IF (i == j) cycle jloop
3247 xij(:) = xpj(:) - xpi(:)
3249 IF (rij > rcut) cycle jloop
3252 CALL accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3253 denominator_tolerance, l, mm, nbond, re_qlm, im_qlm, &
3254 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3264 IF (nbond < denominator_tolerance)
THEN
3265 cpwarn(
"QPARM: number of neighbors is very close to zero!")
3268 d_nbond_dxi(:) = d_nbond_dxi(:)/nbond
3269 re_qlm = re_qlm/nbond
3270 d_re_qlm_dxi(:) = d_re_qlm_dxi(:)/nbond - d_nbond_dxi(:)*re_qlm
3271 im_qlm = im_qlm/nbond
3272 d_im_qlm_dxi(:) = d_im_qlm_dxi(:)/nbond - d_nbond_dxi(:)*im_qlm
3274 ql = ql + fact*(re_qlm*re_qlm + im_qlm*im_qlm)
3275 d_ql_dxi(:) = d_ql_dxi(:) &
3276 + fact*2.0_dp*(re_qlm*d_re_qlm_dxi(:) + im_qlm*d_im_qlm_dxi(:))
3280 pre_fac = (4.0_dp*
pi)/(2.0_dp*l + 1)
3281 qparm = qparm + sqrt(pre_fac*ql)
3282 ftmp(:) = 0.5_dp*sqrt(pre_fac/ql)*d_ql_dxi(:)
3284 ftmp(:) = -1.0_dp*ftmp(:)
3286 CALL put_derivative(colvar, ii, ftmp)
3290 colvar%ss = qparm*inv_n_atoms_from
3291 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
3293 END SUBROUTINE qparm_colvar
3311 SUBROUTINE accumulate_qlm_over_neigbors(xij, rij, rcut, r1cut, &
3312 denominator_tolerance, ll, mm, nbond, re_qlm, im_qlm, &
3313 d_re_qlm_dxi, d_im_qlm_dxi, d_nbond_dxi)
3315 REAL(kind=
dp),
INTENT(IN) :: xij(3), rij, rcut, r1cut, &
3316 denominator_tolerance
3317 INTEGER,
INTENT(IN) :: ll, mm
3318 REAL(kind=
dp),
INTENT(INOUT) :: nbond, re_qlm, im_qlm, d_re_qlm_dxi(3), &
3319 d_im_qlm_dxi(3), d_nbond_dxi(3)
3321 REAL(kind=
dp) :: bond, costheta, dplm, dylm, exp0, &
3322 exp_fac, fi, plm, pre_fac, sqrt_c1
3323 REAL(kind=
dp),
DIMENSION(3) :: dcostheta, dfi
3327 IF (rij > rcut)
THEN
3332 IF (rij < r1cut)
THEN
3336 exp0 = exp((r1cut - rcut)/(rij - rcut) - (r1cut - rcut)/(r1cut - rij))
3337 bond = 1.0_dp/(1.0_dp + exp0)
3338 exp_fac = ((rcut - r1cut)/(rij - rcut)**2 + (rcut - r1cut)/(r1cut - rij)**2)*exp0/(1.0_dp + exp0)**2
3341 IF (bond > 1.0_dp)
THEN
3342 cpabort(
"bond > 1.0_dp")
3345 nbond = nbond + bond
3346 IF (abs(xij(1)) < denominator_tolerance &
3347 .AND. abs(xij(2)) < denominator_tolerance)
THEN
3350 fi = atan2(xij(2), xij(1))
3353 costheta = xij(3)/rij
3354 IF (costheta > 1.0_dp) costheta = 1.0_dp
3355 IF (costheta < -1.0_dp) costheta = -1.0_dp
3360 IF ((ll + abs(mm)) >
maxfac)
THEN
3361 cpabort(
"(l+m) > maxfac")
3364 sqrt_c1 = sqrt(((2*ll + 1)*
fac(ll - abs(mm)))/(4*
pi*
fac(ll + abs(mm))))
3365 pre_fac = bond*sqrt_c1
3368 re_qlm = re_qlm + pre_fac*plm*cos(mm*fi)
3369 im_qlm = im_qlm + pre_fac*plm*sin(mm*fi)
3371 dcostheta(:) = xij(:)*xij(3)/(rij**3)
3372 dcostheta(3) = dcostheta(3) - 1.0_dp/rij
3376 dfi(1) = xij(2)/(xij(1)**2 + xij(2)**2)
3377 dfi(2) = -xij(1)/(xij(1)**2 + xij(2)**2)
3379 d_re_qlm_dxi(:) = d_re_qlm_dxi(:) &
3380 + exp_fac*sqrt_c1*plm*cos(mm*fi)*xij(:)/rij &
3381 + dylm*dcostheta(:)*cos(mm*fi) &
3382 + pre_fac*plm*mm*(-1.0_dp)*sin(mm*fi)*dfi(:)
3383 d_im_qlm_dxi(:) = d_im_qlm_dxi(:) &
3384 + exp_fac*sqrt_c1*plm*sin(mm*fi)*xij(:)/rij &
3385 + dylm*dcostheta(:)*sin(mm*fi) &
3386 + pre_fac*plm*mm*(+1.0_dp)*cos(mm*fi)*dfi(:)
3387 d_nbond_dxi(:) = d_nbond_dxi(:) + exp_fac*xij(:)/rij
3389 END SUBROUTINE accumulate_qlm_over_neigbors
3401 SUBROUTINE hydronium_shell_colvar(colvar, cell, subsys, particles)
3406 POINTER :: particles
3408 INTEGER :: i, ii, j, jj, n_hydrogens, n_oxygens, &
3409 pm, poh, poo, qm, qoh, qoo
3410 REAL(
dp) :: drji, fscalar, invden, lambda, nh, num, &
3411 qtot, rji(3), roh, roo, rrel
3412 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: m, noh, noo, qloc
3413 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dm, dnoh, dnoo
3414 REAL(
dp),
DIMENSION(3) :: rpi, rpj
3418 n_oxygens = colvar%hydronium_shell_param%n_oxygens
3419 n_hydrogens = colvar%hydronium_shell_param%n_hydrogens
3420 nh = colvar%hydronium_shell_param%nh
3421 poh = colvar%hydronium_shell_param%poh
3422 qoh = colvar%hydronium_shell_param%qoh
3423 poo = colvar%hydronium_shell_param%poo
3424 qoo = colvar%hydronium_shell_param%qoo
3425 roo = colvar%hydronium_shell_param%roo
3426 roh = colvar%hydronium_shell_param%roh
3427 lambda = colvar%hydronium_shell_param%lambda
3428 pm = colvar%hydronium_shell_param%pm
3429 qm = colvar%hydronium_shell_param%qm
3431 NULLIFY (particles_i)
3433 IF (
PRESENT(particles))
THEN
3434 my_particles => particles
3436 cpassert(
PRESENT(subsys))
3438 my_particles => particles_i%els
3441 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3442 ALLOCATE (noh(n_oxygens))
3443 ALLOCATE (m(n_oxygens))
3444 ALLOCATE (dm(3, n_hydrogens, n_oxygens))
3446 ALLOCATE (dnoo(3, n_oxygens, n_oxygens))
3447 ALLOCATE (noo(n_oxygens))
3449 ALLOCATE (qloc(n_oxygens))
3459 DO ii = 1, n_oxygens
3460 i = colvar%hydronium_shell_param%i_oxygens(ii)
3461 rpi(:) = my_particles(i)%r(1:3)
3463 DO jj = 1, n_hydrogens
3464 j = colvar%hydronium_shell_param%i_hydrogens(jj)
3465 rpj(:) = my_particles(j)%r(1:3)
3466 rji =
pbc(rpj, rpi, cell)
3467 drji = sqrt(sum(rji**2))
3469 num = (1.0_dp - rrel**poh)
3470 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3471 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3472 noh(ii) = noh(ii) + num*invden
3473 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3474 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3475 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3478 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3479 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3480 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3483 m(ii) = 1.0_dp - (1.0_dp - (noh(ii)/nh)**pm)/ &
3484 (1.0_dp - (noh(ii)/nh)**qm)
3487 DO jj = 1, n_oxygens
3489 j = colvar%hydronium_shell_param%i_oxygens(jj)
3490 rpj(:) = my_particles(j)%r(1:3)
3491 rji =
pbc(rpj, rpi, cell)
3492 drji = sqrt(sum(rji**2))
3494 num = (1.0_dp - rrel**poo)
3495 invden = 1.0_dp/(1.0_dp - rrel**qoo)
3496 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3497 noo(ii) = noo(ii) + num*invden
3498 fscalar = ((-poo*(rrel**(poo - 1))*invden) &
3499 + num*(invden)**2*qoo*(rrel**(qoo - 1)))/(drji*roo)
3500 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3503 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
3504 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drji)
3505 dnoo(1:3, jj, ii) = rji(1:3)*fscalar
3512 DO ii = 1, n_oxygens
3513 qloc(ii) = exp(lambda*m(ii)*noo(ii))
3514 qtot = qtot + qloc(ii)
3517 DO ii = 1, n_oxygens
3519 DO jj = 1, n_hydrogens
3520 dm(1:3, jj, ii) = (pm*((noh(ii)/nh)**(pm - 1))*dnoh(1:3, jj, ii))/nh/ &
3521 (1.0_dp - (noh(ii)/nh)**qm) - &
3522 (1.0_dp - (noh(ii)/nh)**pm)/ &
3523 ((1.0_dp - (noh(ii)/nh)**qm)**2)* &
3524 qm*dnoh(1:3, jj, ii)*(noh(ii)/nh)**(qm - 1)/nh
3526 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3527 colvar%dsdr(1:3, n_oxygens + jj) = colvar%dsdr(1:3, n_oxygens + jj) &
3528 - qloc(ii)*dm(1:3, jj, ii)*noo(ii)/qtot
3531 DO jj = 1, n_oxygens
3532 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3533 colvar%dsdr(1:3, jj) = colvar%dsdr(1:3, jj) &
3534 - qloc(ii)*m(ii)*dnoo(1:3, jj, ii)/qtot
3538 colvar%ss = log(qtot)/lambda
3547 END SUBROUTINE hydronium_shell_colvar
3558 SUBROUTINE hydronium_dist_colvar(colvar, cell, subsys, particles)
3563 POINTER :: particles
3565 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3566 n_oxygens, offseth, pf, pm, poh, qf, &
3568 REAL(
dp) :: drji, drki, fscalar, invden, lambda, nh, nn, num, rion, rion_den, rion_num, &
3569 rji(3), rki(3), roh, rrel, sum_expfac_f, sum_expfac_noh
3570 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac_f, dexpfac_noh, df, dm, &
3571 expfac_f, expfac_f_rki, expfac_noh, f, &
3573 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_f_rki
3574 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rki, dnoh
3575 REAL(
dp),
DIMENSION(3) :: rpi, rpj, rpk
3579 n_oxygens = colvar%hydronium_dist_param%n_oxygens
3580 n_hydrogens = colvar%hydronium_dist_param%n_hydrogens
3581 poh = colvar%hydronium_dist_param%poh
3582 qoh = colvar%hydronium_dist_param%qoh
3583 roh = colvar%hydronium_dist_param%roh
3584 pm = colvar%hydronium_dist_param%pm
3585 qm = colvar%hydronium_dist_param%qm
3586 nh = colvar%hydronium_dist_param%nh
3587 pf = colvar%hydronium_dist_param%pf
3588 qf = colvar%hydronium_dist_param%qf
3589 nn = colvar%hydronium_dist_param%nn
3590 lambda = colvar%hydronium_dist_param%lambda
3592 NULLIFY (particles_i)
3594 IF (
PRESENT(particles))
THEN
3595 my_particles => particles
3597 cpassert(
PRESENT(subsys))
3599 my_particles => particles_i%els
3602 ALLOCATE (dnoh(3, n_hydrogens, n_oxygens))
3603 ALLOCATE (noh(n_oxygens))
3604 ALLOCATE (m(n_oxygens), dm(n_oxygens))
3605 ALLOCATE (f(n_oxygens), df(n_oxygens))
3606 ALLOCATE (expfac_noh(n_oxygens), dexpfac_noh(n_oxygens))
3607 ALLOCATE (expfac_f(n_oxygens), dexpfac_f(n_oxygens))
3608 ALLOCATE (ddist_rki(3, n_oxygens, n_oxygens))
3609 ALLOCATE (expfac_f_rki(n_oxygens))
3610 ALLOCATE (dexpfac_f_rki(n_oxygens, n_oxygens))
3622 sum_expfac_noh = 0._dp
3623 sum_expfac_f = 0._dp
3625 expfac_f_rki = 0._dp
3626 dexpfac_f_rki = 0._dp
3629 DO ii = 1, n_oxygens
3630 i = colvar%hydronium_dist_param%i_oxygens(ii)
3631 rpi(:) = my_particles(i)%r(1:3)
3632 DO jj = 1, n_hydrogens
3633 j = colvar%hydronium_dist_param%i_hydrogens(jj)
3634 rpj(:) = my_particles(j)%r(1:3)
3635 rji =
pbc(rpj, rpi, cell)
3636 drji = sqrt(sum(rji**2))
3638 num = (1.0_dp - rrel**poh)
3639 invden = 1.0_dp/(1.0_dp - rrel**qoh)
3640 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3641 noh(ii) = noh(ii) + num*invden
3642 fscalar = ((-poh*(rrel**(poh - 1))*invden) &
3643 + num*(invden)**2*qoh*(rrel**(qoh - 1)))/(drji*roh)
3644 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3647 noh(ii) = noh(ii) + real(poh,
dp)/real(qoh,
dp)
3648 fscalar = real(poh*(poh - qoh),
dp)/(real(2*qoh,
dp)*roh*drji)
3649 dnoh(1:3, jj, ii) = rji(1:3)*fscalar
3655 DO ii = 1, n_oxygens
3656 num = 1.0_dp - (noh(ii)/nh)**pm
3657 invden = 1.0_dp/(1.0_dp - (noh(ii)/nh)**qm)
3658 m(ii) = 1.0_dp - num*invden
3659 dm(ii) = (pm*(noh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
3660 (noh(ii)/nh)**(qm - 1))/nh
3661 expfac_noh(ii) = exp(lambda*noh(ii))
3662 dexpfac_noh(ii) = lambda*expfac_noh(ii)
3663 sum_expfac_noh = sum_expfac_noh + expfac_noh(ii)
3667 DO ii = 1, n_oxygens
3668 i = colvar%hydronium_dist_param%i_oxygens(ii)
3669 num = 1.0_dp - (noh(ii)/nn)**pf
3670 invden = 1.0_dp/(1.0_dp - (noh(ii)/nn)**qf)
3672 df(ii) = (-pf*(noh(ii)/nn)**(pf - 1)*invden + qf*num*(invden**2)* &
3673 (noh(ii)/nn)**(qf - 1))/nn
3674 expfac_f(ii) = exp(lambda*f(ii))
3675 dexpfac_f(ii) = lambda*expfac_f(ii)
3676 sum_expfac_f = sum_expfac_f + expfac_f(ii)
3680 DO ii = 1, n_oxygens
3681 i = colvar%hydronium_dist_param%i_oxygens(ii)
3682 rpi(:) = my_particles(i)%r(1:3)
3683 DO kk = 1, n_oxygens
3685 k = colvar%hydronium_dist_param%i_oxygens(kk)
3686 rpk(:) = my_particles(k)%r(1:3)
3687 rki =
pbc(rpk, rpi, cell)
3688 drki = sqrt(sum(rki**2))
3689 expfac_f_rki(ii) = expfac_f_rki(ii) + drki*expfac_f(kk)
3690 ddist_rki(1:3, kk, ii) = rki(1:3)/drki
3691 dexpfac_f_rki(kk, ii) = drki*dexpfac_f(kk)
3693 rion_num = rion_num + m(ii)*expfac_noh(ii)*expfac_f_rki(ii)
3697 rion_den = sum_expfac_noh*sum_expfac_f
3698 rion = rion_num/rion_den
3703 DO ii = 1, n_oxygens
3704 DO jj = 1, n_hydrogens
3705 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3706 + dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3707 *expfac_f_rki(ii)/rion_den
3708 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3709 - dm(ii)*dnoh(1:3, jj, ii)*expfac_noh(ii) &
3710 *expfac_f_rki(ii)/rion_den
3711 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3712 + m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3713 *expfac_f_rki(ii)/rion_den
3714 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3715 - m(ii)*dexpfac_noh(ii)*dnoh(1:3, jj, ii) &
3716 *expfac_f_rki(ii)/rion_den
3718 DO kk = 1, n_oxygens
3720 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3721 - m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3722 *expfac_f(kk)/rion_den
3723 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3724 + m(ii)*expfac_noh(ii)*ddist_rki(1:3, kk, ii) &
3725 *expfac_f(kk)/rion_den
3726 DO jj = 1, n_hydrogens
3727 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) &
3728 + m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3729 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3730 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3731 - m(ii)*expfac_noh(ii)*dexpfac_f_rki(kk, ii) &
3732 *df(kk)*dnoh(1:3, jj, kk)/rion_den
3737 DO ii = 1, n_oxygens
3738 DO jj = 1, n_hydrogens
3739 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3740 - rion_num*sum_expfac_f*dexpfac_noh(ii) &
3741 *dnoh(1:3, jj, ii)/(rion_den**2)
3742 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3743 + rion_num*sum_expfac_f*dexpfac_noh(ii) &
3744 *dnoh(1:3, jj, ii)/(rion_den**2)
3745 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3746 - rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3747 *dnoh(1:3, jj, ii)/(rion_den**2)
3748 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3749 + rion_num*sum_expfac_noh*dexpfac_f(ii)*df(ii) &
3750 *dnoh(1:3, jj, ii)/(rion_den**2)
3754 DEALLOCATE (noh, m, f, expfac_noh, expfac_f)
3755 DEALLOCATE (dnoh, dm, df, dexpfac_noh, dexpfac_f)
3756 DEALLOCATE (ddist_rki, expfac_f_rki, dexpfac_f_rki)
3758 END SUBROUTINE hydronium_dist_colvar
3771 SUBROUTINE acid_hyd_dist_colvar(colvar, cell, subsys, particles)
3776 POINTER :: particles
3778 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, &
3779 n_oxygens_acid, n_oxygens_water, &
3780 offseth, offseto, paoh, pcut, pwoh, &
3782 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dexpfac, expfac, nwoh
3783 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dexpfac_rik
3784 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ddist_rik, dnaoh, dnwoh
3785 REAL(kind=
dp) :: dfcut, drik, drji, drjk, fbrace, fcut, fscalar, invden, invden_cut, lambda, &
3786 naoh, nc, num, num_cut, raoh, rik(3), rion, rion_den, rion_num, rji(3), rjk(3), rpi(3), &
3787 rpj(3), rpk(3), rrel, rwoh
3791 NULLIFY (my_particles, particles_i)
3793 n_oxygens_water = colvar%acid_hyd_dist_param%n_oxygens_water
3794 n_oxygens_acid = colvar%acid_hyd_dist_param%n_oxygens_acid
3795 n_hydrogens = colvar%acid_hyd_dist_param%n_hydrogens
3796 pwoh = colvar%acid_hyd_dist_param%pwoh
3797 qwoh = colvar%acid_hyd_dist_param%qwoh
3798 paoh = colvar%acid_hyd_dist_param%paoh
3799 qaoh = colvar%acid_hyd_dist_param%qaoh
3800 pcut = colvar%acid_hyd_dist_param%pcut
3801 qcut = colvar%acid_hyd_dist_param%qcut
3802 rwoh = colvar%acid_hyd_dist_param%rwoh
3803 raoh = colvar%acid_hyd_dist_param%raoh
3804 nc = colvar%acid_hyd_dist_param%nc
3805 lambda = colvar%acid_hyd_dist_param%lambda
3806 ALLOCATE (expfac(n_oxygens_water))
3807 ALLOCATE (nwoh(n_oxygens_water))
3808 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
3809 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
3810 ALLOCATE (dexpfac(n_oxygens_water))
3811 ALLOCATE (ddist_rik(3, n_oxygens_water, n_oxygens_acid))
3812 ALLOCATE (dexpfac_rik(n_oxygens_water, n_oxygens_acid))
3817 dnaoh(:, :, :) = 0._dp
3818 dnwoh(:, :, :) = 0._dp
3819 ddist_rik(:, :, :) = 0._dp
3821 dexpfac_rik(:, :) = 0._dp
3824 IF (
PRESENT(particles))
THEN
3825 my_particles => particles
3827 cpassert(
PRESENT(subsys))
3829 my_particles => particles_i%els
3833 DO ii = 1, n_oxygens_water
3834 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3835 rpi(:) = my_particles(i)%r(1:3)
3836 DO jj = 1, n_hydrogens
3837 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3838 rpj(:) = my_particles(j)%r(1:3)
3839 rji =
pbc(rpj, rpi, cell)
3840 drji = sqrt(sum(rji**2))
3842 num = 1.0_dp - rrel**pwoh
3843 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
3844 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3845 nwoh(ii) = nwoh(ii) + num*invden
3846 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
3847 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
3848 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3851 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
3852 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
3853 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
3856 expfac(ii) = exp(lambda*nwoh(ii))
3857 dexpfac(ii) = lambda*expfac(ii)
3858 rion_den = rion_den + expfac(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 ii = 1, n_oxygens_water
3866 i = colvar%acid_hyd_dist_param%i_oxygens_water(ii)
3867 rpi(:) = my_particles(i)%r(1:3)
3868 rik =
pbc(rpi, rpk, cell)
3869 drik = sqrt(sum(rik**2))
3870 rion_num = rion_num + drik*expfac(ii)
3871 ddist_rik(1:3, ii, kk) = rik(1:3)/drik
3872 dexpfac_rik(ii, kk) = drik*dexpfac(ii)
3877 DO kk = 1, n_oxygens_acid
3878 k = colvar%acid_hyd_dist_param%i_oxygens_acid(kk)
3879 rpk(:) = my_particles(k)%r(1:3)
3880 DO jj = 1, n_hydrogens
3881 j = colvar%acid_hyd_dist_param%i_hydrogens(jj)
3882 rpj(:) = my_particles(j)%r(1:3)
3883 rjk =
pbc(rpj, rpk, cell)
3884 drjk = sqrt(sum(rjk**2))
3886 num = 1.0_dp - rrel**paoh
3887 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
3888 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
3889 naoh = naoh + num*invden
3890 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
3891 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
3892 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3895 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
3896 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
3897 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
3901 num_cut = 1.0_dp - (naoh/nc)**pcut
3902 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
3903 fcut = num_cut*invden_cut
3907 fbrace = rion_num/rion_den/n_oxygens_acid
3912 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
3913 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
3914 offseto = n_oxygens_water
3915 offseth = n_oxygens_water + n_oxygens_acid
3916 DO kk = 1, n_oxygens_acid
3917 DO jj = 1, n_hydrogens
3918 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3919 + dfcut*dnaoh(1:3, jj, kk)*fbrace
3920 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3921 - dfcut*dnaoh(1:3, jj, kk)*fbrace
3927 DO kk = 1, n_oxygens_acid
3928 DO ii = 1, n_oxygens_water
3929 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
3930 + fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3932 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3933 - fcut*ddist_rik(1:3, ii, kk)*expfac(ii)/rion_den/n_oxygens_acid
3935 DO jj = 1, n_hydrogens
3936 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3937 + fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3939 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3940 - fcut*dexpfac_rik(ii, kk)*dnwoh(1:3, jj, ii)/rion_den/n_oxygens_acid
3946 DO ii = 1, n_oxygens_water
3947 DO jj = 1, n_hydrogens
3948 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
3949 - fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3950 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
3951 + fcut*rion_num*dexpfac(ii)*dnwoh(1:3, jj, ii)/2.0_dp/(rion_den**2)
3955 END SUBROUTINE acid_hyd_dist_colvar
3968 SUBROUTINE acid_hyd_shell_colvar(colvar, cell, subsys, particles)
3973 POINTER :: particles
3975 INTEGER :: i, ii, j, jj, k, kk, n_hydrogens, n_oxygens_acid, n_oxygens_water, offseth, &
3976 offseto, paoh, pcut, pm, poo, pwoh, qaoh, qcut, qm, qoo, qwoh, tt
3977 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dm, m, noo, nwoh, qloc
3978 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dnaoh, dnoo, dnwoh
3979 REAL(kind=
dp) :: dfcut, drji, drjk, drki, fcut, fscalar, invden, invden_cut, lambda, naoh, &
3980 nc, nh, num, num_cut, qsol, qtot, raoh, rji(3), rjk(3), rki(3), roo, rpi(3), rpj(3), &
3985 NULLIFY (my_particles, particles_i)
3987 n_oxygens_water = colvar%acid_hyd_shell_param%n_oxygens_water
3988 n_oxygens_acid = colvar%acid_hyd_shell_param%n_oxygens_acid
3989 n_hydrogens = colvar%acid_hyd_shell_param%n_hydrogens
3990 pwoh = colvar%acid_hyd_shell_param%pwoh
3991 qwoh = colvar%acid_hyd_shell_param%qwoh
3992 paoh = colvar%acid_hyd_shell_param%paoh
3993 qaoh = colvar%acid_hyd_shell_param%qaoh
3994 poo = colvar%acid_hyd_shell_param%poo
3995 qoo = colvar%acid_hyd_shell_param%qoo
3996 pm = colvar%acid_hyd_shell_param%pm
3997 qm = colvar%acid_hyd_shell_param%qm
3998 pcut = colvar%acid_hyd_shell_param%pcut
3999 qcut = colvar%acid_hyd_shell_param%qcut
4000 rwoh = colvar%acid_hyd_shell_param%rwoh
4001 raoh = colvar%acid_hyd_shell_param%raoh
4002 roo = colvar%acid_hyd_shell_param%roo
4003 nc = colvar%acid_hyd_shell_param%nc
4004 nh = colvar%acid_hyd_shell_param%nh
4005 lambda = colvar%acid_hyd_shell_param%lambda
4006 ALLOCATE (nwoh(n_oxygens_water))
4007 ALLOCATE (dnwoh(3, n_hydrogens, n_oxygens_water))
4008 ALLOCATE (dnaoh(3, n_hydrogens, n_oxygens_acid))
4009 ALLOCATE (m(n_oxygens_water))
4010 ALLOCATE (dm(n_oxygens_water))
4011 ALLOCATE (noo(n_oxygens_water))
4012 ALLOCATE (dnoo(3, n_oxygens_water + n_oxygens_acid, n_oxygens_water))
4013 ALLOCATE (qloc(n_oxygens_water))
4017 dnaoh(:, :, :) = 0._dp
4018 dnwoh(:, :, :) = 0._dp
4019 dnoo(:, :, :) = 0._dp
4025 IF (
PRESENT(particles))
THEN
4026 my_particles => particles
4028 cpassert(
PRESENT(subsys))
4030 my_particles => particles_i%els
4034 DO ii = 1, n_oxygens_water
4035 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4036 rpi(:) = my_particles(i)%r(1:3)
4037 DO jj = 1, n_hydrogens
4038 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4039 rpj(:) = my_particles(j)%r(1:3)
4040 rji =
pbc(rpj, rpi, cell)
4041 drji = sqrt(sum(rji**2))
4043 num = 1.0_dp - rrel**pwoh
4044 invden = 1.0_dp/(1.0_dp - rrel**qwoh)
4045 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4046 nwoh(ii) = nwoh(ii) + num*invden
4047 fscalar = (-pwoh*(rrel**(pwoh - 1))*invden &
4048 + num*(invden**2)*qwoh*(rrel**(qwoh - 1)))/(drji*rwoh)
4049 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4052 nwoh(ii) = nwoh(ii) + real(pwoh,
dp)/real(qwoh,
dp)
4053 fscalar = real(pwoh*(pwoh - qwoh),
dp)/(real(2*qwoh,
dp)*rwoh*drji)
4054 dnwoh(1:3, jj, ii) = rji(1:3)*fscalar
4060 DO ii = 1, n_oxygens_water
4061 num = 1.0_dp - (nwoh(ii)/nh)**pm
4062 invden = 1.0_dp/(1.0_dp - (nwoh(ii)/nh)**qm)
4063 m(ii) = 1.0_dp - num*invden
4064 dm(ii) = (pm*(nwoh(ii)/nh)**(pm - 1)*invden - qm*num*(invden**2)* &
4065 (nwoh(ii)/nh)**(qm - 1))/nh
4069 DO ii = 1, n_oxygens_water
4070 i = colvar%acid_hyd_shell_param%i_oxygens_water(ii)
4071 rpi(:) = my_particles(i)%r(1:3)
4072 DO kk = 1, n_oxygens_water + n_oxygens_acid
4074 IF (kk <= n_oxygens_water)
THEN
4075 k = colvar%acid_hyd_shell_param%i_oxygens_water(kk)
4076 rpk(:) = my_particles(k)%r(1:3)
4078 tt = kk - n_oxygens_water
4079 k = colvar%acid_hyd_shell_param%i_oxygens_acid(tt)
4080 rpk(:) = my_particles(k)%r(1:3)
4082 rki =
pbc(rpk, rpi, cell)
4083 drki = sqrt(sum(rki**2))
4085 num = 1.0_dp - rrel**poo
4086 invden = 1.0_dp/(1.0_dp - rrel**qoo)
4087 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4088 noo(ii) = noo(ii) + num*invden
4089 fscalar = (-poo*(rrel**(poo - 1))*invden &
4090 + num*(invden**2)*qoo*(rrel**(qoo - 1)))/(drki*roo)
4091 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4094 noo(ii) = noo(ii) + real(poo,
dp)/real(qoo,
dp)
4095 fscalar = real(poo*(poo - qoo),
dp)/(real(2*qoo,
dp)*roo*drki)
4096 dnoo(1:3, kk, ii) = rki(1:3)*fscalar
4102 DO kk = 1, n_oxygens_acid
4103 k = colvar%acid_hyd_shell_param%i_oxygens_acid(kk)
4104 rpk(:) = my_particles(k)%r(1:3)
4105 DO jj = 1, n_hydrogens
4106 j = colvar%acid_hyd_shell_param%i_hydrogens(jj)
4107 rpj(:) = my_particles(j)%r(1:3)
4108 rjk =
pbc(rpj, rpk, cell)
4109 drjk = sqrt(sum(rjk**2))
4111 num = 1.0_dp - rrel**paoh
4112 invden = 1.0_dp/(1.0_dp - rrel**qaoh)
4113 IF (abs(1.0_dp - rrel) > 1.0e-6_dp)
THEN
4114 naoh = naoh + num*invden
4115 fscalar = (-paoh*(rrel**(paoh - 1))*invden &
4116 + num*(invden**2)*qaoh*(rrel**(qaoh - 1)))/(drjk*raoh)
4117 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4120 naoh = naoh + real(paoh,
dp)/real(qaoh,
dp)
4121 fscalar = real(paoh*(paoh - qaoh),
dp)/(real(2*qaoh,
dp)*raoh*drjk)
4122 dnaoh(1:3, jj, kk) = rjk(1:3)*fscalar
4126 num_cut = 1.0_dp - (naoh/nc)**pcut
4127 invden_cut = 1.0_dp/(1.0_dp - (naoh/nc)**qcut)
4128 fcut = num_cut*invden_cut
4131 DO ii = 1, n_oxygens_water
4132 qloc(ii) = exp(lambda*m(ii)*noo(ii))
4133 qtot = qtot + qloc(ii)
4135 qsol = log(qtot)/lambda
4136 colvar%ss = fcut*qsol
4139 dfcut = ((-pcut*(naoh/nc)**(pcut - 1)*invden_cut) &
4140 + num_cut*(invden_cut**2)*qcut*(naoh/nc)**(qcut - 1))/nc
4141 offseto = n_oxygens_water
4142 offseth = n_oxygens_water + n_oxygens_acid
4143 DO kk = 1, n_oxygens_acid
4144 DO jj = 1, n_hydrogens
4145 colvar%dsdr(1:3, offseto + kk) = colvar%dsdr(1:3, offseto + kk) &
4146 + dfcut*dnaoh(1:3, jj, kk)*qsol
4147 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4148 - dfcut*dnaoh(1:3, jj, kk)*qsol
4154 DO ii = 1, n_oxygens_water
4155 fscalar = fcut*qloc(ii)*dm(ii)*noo(ii)/qtot
4156 DO jj = 1, n_hydrogens
4157 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) &
4158 + fscalar*dnwoh(1:3, jj, ii)
4159 colvar%dsdr(1:3, offseth + jj) = colvar%dsdr(1:3, offseth + jj) &
4160 - fscalar*dnwoh(1:3, jj, ii)
4164 DO ii = 1, n_oxygens_water
4165 fscalar = fcut*qloc(ii)*m(ii)/qtot
4166 DO kk = 1, n_oxygens_water + n_oxygens_acid
4168 colvar%dsdr(1:3, ii) = colvar%dsdr(1:3, ii) + fscalar*dnoo(1:3, kk, ii)
4169 colvar%dsdr(1:3, kk) = colvar%dsdr(1:3, kk) - fscalar*dnoo(1:3, kk, ii)
4173 END SUBROUTINE acid_hyd_shell_colvar
4185 SUBROUTINE coord_colvar(colvar, cell, subsys, particles)
4190 POINTER :: particles
4192 INTEGER :: i, ii, j, jj, k, kk, n_atoms_from, &
4193 n_atoms_to_a, n_atoms_to_b, p_a, p_b, &
4195 REAL(
dp) :: dfunc_ij, dfunc_jk, func_ij, func_jk, func_k, inv_n_atoms_from, invden_ij, &
4196 invden_jk, ncoord, num_ij, num_jk, r_0_a, r_0_b, rdist_ij, rdist_jk, rij, rjk
4197 REAL(
dp),
DIMENSION(3) :: ftmp_i, ftmp_j, ftmp_k, ss, xij, xjk, &
4205 NULLIFY (particles_i)
4207 IF (
PRESENT(particles))
THEN
4208 my_particles => particles
4210 cpassert(
PRESENT(subsys))
4212 my_particles => particles_i%els
4214 n_atoms_to_a = colvar%coord_param%n_atoms_to
4215 n_atoms_to_b = colvar%coord_param%n_atoms_to_b
4216 n_atoms_from = colvar%coord_param%n_atoms_from
4217 p_a = colvar%coord_param%nncrd
4218 q_a = colvar%coord_param%ndcrd
4219 r_0_a = colvar%coord_param%r_0
4220 p_b = colvar%coord_param%nncrd_b
4221 q_b = colvar%coord_param%ndcrd_b
4222 r_0_b = colvar%coord_param%r_0_b
4225 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
4226 DO ii = 1, n_atoms_from
4227 i = colvar%coord_param%i_at_from(ii)
4228 CALL get_coordinates(colvar, i, xpi, my_particles)
4229 DO jj = 1, n_atoms_to_a
4230 j = colvar%coord_param%i_at_to(jj)
4231 CALL get_coordinates(colvar, j, xpj, my_particles)
4234 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
4236 xij = matmul(cell%hmat, ss)
4237 rij = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
4238 IF (rij < 1.0e-8_dp) cycle
4239 rdist_ij = rij/r_0_a
4240 IF (abs(1.0_dp - rdist_ij) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4241 num_ij = (1.0_dp - rdist_ij**p_a)
4242 invden_ij = 1.0_dp/(1.0_dp - rdist_ij**q_a)
4243 func_ij = num_ij*invden_ij
4244 IF (rij < 1.0e-8_dp)
THEN
4248 dfunc_ij = (-p_a*rdist_ij**(p_a - 1)*invden_ij &
4249 + num_ij*(invden_ij)**2*q_a*rdist_ij**(q_a - 1))/(rij*r_0_a)
4253 func_ij = real(p_a, kind=
dp)/real(q_a, kind=
dp)
4254 dfunc_ij = real(p_a, kind=
dp)*real((-q_a + p_a), kind=
dp)/(real(2*q_a, kind=
dp)*r_0_a)
4256 IF (n_atoms_to_b /= 0)
THEN
4258 DO kk = 1, n_atoms_to_b
4259 k = colvar%coord_param%i_at_to_b(kk)
4261 CALL get_coordinates(colvar, k, xpk, my_particles)
4262 ss = matmul(cell%h_inv, xpj(:) - xpk(:))
4264 xjk = matmul(cell%hmat, ss)
4265 rjk = sqrt(xjk(1)**2 + xjk(2)**2 + xjk(3)**2)
4266 IF (rjk < 1.0e-8_dp) cycle
4267 rdist_jk = rjk/r_0_b
4268 IF (abs(1.0_dp - rdist_jk) > epsilon(0.0_dp)*1.0e+4_dp)
THEN
4269 num_jk = (1.0_dp - rdist_jk**p_b)
4270 invden_jk = 1.0_dp/(1.0_dp - rdist_jk**q_b)
4271 func_jk = num_jk*invden_jk
4272 IF (rjk < 1.0e-8_dp)
THEN
4276 dfunc_jk = (-p_b*rdist_jk**(p_b - 1)*invden_jk &
4277 + num_jk*(invden_jk)**2*q_b*rdist_jk**(q_b - 1))/(rjk*r_0_b)
4281 func_jk = real(p_b, kind=
dp)/real(q_b, kind=
dp)
4282 dfunc_jk = real(p_b, kind=
dp)*real((-q_b + p_b), kind=
dp)/(real(2*q_b, kind=
dp)*r_0_b)
4284 func_k = func_k + func_jk
4285 ftmp_k = -func_ij*dfunc_jk*xjk
4286 CALL put_derivative(colvar, n_atoms_from + n_atoms_to_a + kk, ftmp_k)
4288 ftmp_j = -dfunc_ij*xij*func_jk + func_ij*dfunc_jk*xjk
4289 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4294 ftmp_j = -dfunc_ij*xij
4295 CALL put_derivative(colvar, n_atoms_from + jj, ftmp_j)
4297 ncoord = ncoord + func_ij*func_k
4298 ftmp_i = dfunc_ij*xij*func_k
4299 CALL put_derivative(colvar, ii, ftmp_i)
4302 colvar%ss = ncoord*inv_n_atoms_from
4303 colvar%dsdr(:, :) = colvar%dsdr(:, :)*inv_n_atoms_from
4304 END SUBROUTINE coord_colvar
4313 SUBROUTINE mindist_colvar(colvar, cell, subsys, particles)
4319 POINTER :: particles
4321 INTEGER :: i, ii, j, jj, n_coord_from, n_coord_to, &
4323 REAL(
dp) :: den_n, den_q, fscalar, ftemp_i(3), inv_den_n, inv_den_q, lambda, num_n, num_q, &
4324 qfunc, r12, r_cut, rfact, rij(3), rpi(3), rpj(3)
4325 REAL(
dp),
DIMENSION(:),
POINTER :: dqfunc_dnl, expnl, nlcoord, sum_rij
4326 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: dnlcoord, dqfunc_dr
4333 NULLIFY (particles_i)
4335 IF (
PRESENT(particles))
THEN
4336 my_particles => particles
4338 cpassert(
PRESENT(subsys))
4340 my_particles => particles_i%els
4343 n_dist_from = colvar%mindist_param%n_dist_from
4344 n_coord_from = colvar%mindist_param%n_coord_from
4345 n_coord_to = colvar%mindist_param%n_coord_to
4346 p = colvar%mindist_param%p_exp
4347 q = colvar%mindist_param%q_exp
4348 r_cut = colvar%mindist_param%r_cut
4349 lambda = colvar%mindist_param%lambda
4351 NULLIFY (nlcoord, dnlcoord, dqfunc_dr, dqfunc_dnl, expnl, sum_rij)
4352 ALLOCATE (nlcoord(n_coord_from))
4353 ALLOCATE (dnlcoord(3, n_coord_from, n_coord_to))
4354 ALLOCATE (expnl(n_coord_from))
4355 ALLOCATE (sum_rij(n_coord_from))
4356 ALLOCATE (dqfunc_dr(3, n_dist_from, n_coord_from))
4357 ALLOCATE (dqfunc_dnl(n_coord_from))
4364 DO i = 1, n_coord_from
4365 ii = colvar%mindist_param%i_coord_from(i)
4366 rpi = my_particles(ii)%r(1:3)
4367 DO j = 1, n_coord_to
4368 jj = colvar%mindist_param%i_coord_to(j)
4369 rpj = my_particles(jj)%r(1:3)
4370 rij =
pbc(rpj, rpi, cell)
4371 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4373 num_n = 1.0_dp - rfact**p
4374 den_n = 1.0_dp - rfact**q
4375 inv_den_n = 1.0_dp/den_n
4376 IF (abs(inv_den_n) < 1.e-10_dp)
THEN
4377 inv_den_n = 1.e-10_dp
4381 fscalar = (-p*rfact**(p - 1) + num_n*q*rfact**(q - 1)*inv_den_n)*inv_den_n/(r_cut*r12)
4383 dnlcoord(1, i, j) = rij(1)*fscalar
4384 dnlcoord(2, i, j) = rij(2)*fscalar
4385 dnlcoord(3, i, j) = rij(3)*fscalar
4387 nlcoord(i) = nlcoord(i) + num_n*inv_den_n
4389 expnl(i) = exp(lambda*nlcoord(i))
4390 den_q = den_q + expnl(i)
4392 inv_den_q = 1.0_dp/den_q
4399 DO i = 1, n_dist_from
4400 ii = colvar%mindist_param%i_dist_from(i)
4401 rpi = my_particles(ii)%r(1:3)
4402 DO j = 1, n_coord_from
4403 jj = colvar%mindist_param%i_coord_from(j)
4404 rpj = my_particles(jj)%r(1:3)
4405 rij =
pbc(rpj, rpi, cell)
4406 r12 = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
4408 num_q = num_q + r12*expnl(j)
4410 sum_rij(j) = sum_rij(j) + r12
4411 dqfunc_dr(1, i, j) = expnl(j)*rij(1)/r12
4412 dqfunc_dr(2, i, j) = expnl(j)*rij(2)/r12
4413 dqfunc_dr(3, i, j) = expnl(j)*rij(3)/r12
4420 qfunc = num_q*inv_den_q
4421 dqfunc_dr = dqfunc_dr*inv_den_q
4424 DO i = 1, n_coord_from
4425 dqfunc_dnl(i) = lambda*expnl(i)*inv_den_q*(sum_rij(i) - num_q*inv_den_q)
4429 DO i = 1, n_dist_from
4430 DO j = 1, n_coord_from
4431 ftemp_i(1) = dqfunc_dr(1, i, j)
4432 ftemp_i(2) = dqfunc_dr(2, i, j)
4433 ftemp_i(3) = dqfunc_dr(3, i, j)
4435 CALL put_derivative(colvar, i, ftemp_i)
4436 CALL put_derivative(colvar, j + n_dist_from, -ftemp_i)
4440 DO i = 1, n_coord_from
4441 DO j = 1, n_coord_to
4442 ftemp_i(1) = dqfunc_dnl(i)*dnlcoord(1, i, j)
4443 ftemp_i(2) = dqfunc_dnl(i)*dnlcoord(2, i, j)
4444 ftemp_i(3) = dqfunc_dnl(i)*dnlcoord(3, i, j)
4446 CALL put_derivative(colvar, i + n_dist_from, ftemp_i)
4447 CALL put_derivative(colvar, j + n_dist_from + n_coord_from, -ftemp_i)
4452 DEALLOCATE (nlcoord)
4453 DEALLOCATE (dnlcoord)
4455 DEALLOCATE (dqfunc_dr)
4456 DEALLOCATE (sum_rij)
4457 DEALLOCATE (dqfunc_dnl)
4459 END SUBROUTINE mindist_colvar
4469 SUBROUTINE combine_colvar(colvar, cell, subsys, particles)
4474 POINTER :: particles
4476 CHARACTER(LEN=default_string_length) :: def_error, this_error
4477 CHARACTER(LEN=default_string_length), &
4478 ALLOCATABLE,
DIMENSION(:) :: my_par
4479 INTEGER :: i, ii, j, ncolv, ndim
4481 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: dss_vals, my_val, ss_vals
4482 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi
4487 IF (
PRESENT(particles))
THEN
4488 my_particles => particles
4490 cpassert(
PRESENT(subsys))
4492 my_particles => particles_i%els
4495 ncolv =
SIZE(colvar%combine_cvs_param%colvar_p)
4496 ALLOCATE (ss_vals(ncolv))
4497 ALLOCATE (dss_vals(ncolv))
4501 CALL colvar_recursive_eval(colvar%combine_cvs_param%colvar_p(i)%colvar, cell, my_particles)
4502 ss_vals(i) = colvar%combine_cvs_param%colvar_p(i)%colvar%ss
4507 ndim =
SIZE(colvar%combine_cvs_param%c_parameters) + &
4508 SIZE(colvar%combine_cvs_param%variables)
4509 ALLOCATE (my_par(ndim))
4510 my_par(1:
SIZE(colvar%combine_cvs_param%variables)) = colvar%combine_cvs_param%variables
4511 my_par(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%c_parameters
4512 ALLOCATE (my_val(ndim))
4513 my_val(1:
SIZE(colvar%combine_cvs_param%variables)) = ss_vals
4514 my_val(
SIZE(colvar%combine_cvs_param%variables) + 1:) = colvar%combine_cvs_param%v_parameters
4515 CALL parsef(1, trim(colvar%combine_cvs_param%function), my_par)
4516 colvar%ss =
evalf(1, my_val)
4518 dss_vals(i) =
evalfd(1, i, my_val, colvar%combine_cvs_param%dx, err)
4519 IF ((abs(err) > colvar%combine_cvs_param%lerr))
THEN
4520 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
4521 WRITE (def_error,
"(A,G12.6,A)")
"(", colvar%combine_cvs_param%lerr,
")"
4524 CALL cp_warn(__location__, &
4525 'ASSERTION (cond) failed at line '//
cp_to_string(__line__)// &
4526 ' Error '//trim(this_error)//
' in computing numerical derivatives larger than'// &
4527 trim(def_error)//
' . ')
4535 ALLOCATE (fi(3, colvar%n_atom_s))
4538 DO j = 1, colvar%combine_cvs_param%colvar_p(i)%colvar%n_atom_s
4540 fi(:, ii) = colvar%combine_cvs_param%colvar_p(i)%colvar%dsdr(:, j)*dss_vals(i)
4544 DO i = 1, colvar%n_atom_s
4545 CALL put_derivative(colvar, i, fi(:, i))
4549 DEALLOCATE (ss_vals)
4550 DEALLOCATE (dss_vals)
4551 END SUBROUTINE combine_colvar
4567 SUBROUTINE reaction_path_colvar(colvar, cell, subsys, particles)
4572 POINTER :: particles
4578 IF (
PRESENT(particles))
THEN
4579 my_particles => particles
4581 cpassert(
PRESENT(subsys))
4583 my_particles => particles_i%els
4586 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4587 CALL rpath_dist_rmsd(colvar, my_particles)
4588 ELSE IF (colvar%reaction_path_param%rmsd)
THEN
4589 CALL rpath_rmsd(colvar, my_particles)
4591 CALL rpath_colvar(colvar, cell, my_particles)
4594 END SUBROUTINE reaction_path_colvar
4605 SUBROUTINE rpath_colvar(colvar, cell, particles)
4610 INTEGER :: i, iend, ii, istart, j, k, ncolv, nconf
4611 REAL(
dp) :: lambda, step_size
4612 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: s1, ss_vals
4613 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, f_vals, fi, s1v
4614 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
4616 istart = colvar%reaction_path_param%function_bounds(1)
4617 iend = colvar%reaction_path_param%function_bounds(2)
4619 nconf = colvar%reaction_path_param%nr_frames
4620 step_size = colvar%reaction_path_param%step_size
4621 ncolv = colvar%reaction_path_param%n_components
4622 lambda = colvar%reaction_path_param%lambda
4623 ALLOCATE (f_vals(ncolv, istart:iend))
4624 f_vals(:, :) = colvar%reaction_path_param%f_vals
4625 ALLOCATE (ss_vals(ncolv))
4628 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4629 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4632 ALLOCATE (s1v(2, istart:iend))
4633 ALLOCATE (ds1v(ncolv, 2, istart:iend))
4636 ALLOCATE (ds1(ncolv, 2))
4639 s1v(1, k) = real(k, kind=
dp)*step_size*exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4640 s1v(2, k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4642 ds1v(j, 1, k) = f_vals(j, k)*s1v(1, k)
4643 ds1v(j, 2, k) = f_vals(j, k)*s1v(2, k)
4653 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4655 ALLOCATE (fi(3, colvar%n_atom_s))
4659 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
4661 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)*lambda* &
4662 (ds1(i, 1)/s1(2)/real(nconf - 1,
dp) - colvar%ss*ds1(i, 2)/s1(2))*2.0_dp
4666 DO i = 1, colvar%n_atom_s
4667 CALL put_derivative(colvar, i, fi(:, i))
4672 DEALLOCATE (ss_vals)
4678 END SUBROUTINE rpath_colvar
4689 SUBROUTINE rpath_dist_rmsd(colvar, particles)
4693 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4694 INTEGER,
DIMENSION(:),
POINTER :: iatom
4695 REAL(
dp) :: lambda, my_rmsd, s1(2), sum_exp
4696 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, vec_dif
4697 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dvec_dif, fi, riat, s1v
4698 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4699 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4700 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4702 nconf = colvar%reaction_path_param%nr_frames
4703 rmsd_atom = colvar%reaction_path_param%n_components
4704 lambda = colvar%reaction_path_param%lambda
4705 path_conf => colvar%reaction_path_param%r_ref
4706 iatom => colvar%reaction_path_param%i_rmsd
4708 natom =
SIZE(particles)
4710 ALLOCATE (r0(3*natom))
4711 ALLOCATE (r(3*natom))
4712 ALLOCATE (riat(3, rmsd_atom))
4713 ALLOCATE (vec_dif(rmsd_atom))
4714 ALLOCATE (dvec_dif(3, rmsd_atom))
4715 ALLOCATE (s1v(2, nconf))
4716 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4717 ALLOCATE (ds1(3, rmsd_atom, 2))
4720 r0(ii + 1) = particles(i)%r(1)
4721 r0(ii + 2) = particles(i)%r(2)
4722 r0(ii + 3) = particles(i)%r(3)
4725 DO iat = 1, rmsd_atom
4727 riat(:, iat) = particles(ii)%r
4733 r(ii + 1) = path_conf(ii + 1, ik)
4734 r(ii + 2) = path_conf(ii + 2, ik)
4735 r(ii + 3) = path_conf(ii + 3, ik)
4738 CALL rmsd3(particles, r, r0, output_unit=-1, my_val=my_rmsd, rotate=.true.)
4741 DO iat = 1, rmsd_atom
4744 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 &
4745 + (riat(3, iat) - r(ii + 3))**2
4746 sum_exp = sum_exp + vec_dif(iat)
4749 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*sum_exp)
4750 s1v(2, ik) = exp(-lambda*sum_exp)
4751 DO iat = 1, rmsd_atom
4754 ds1v(1, iat, 1, ik) = r(ii + 1)*s1v(1, ik)
4755 ds1v(1, iat, 2, ik) = r(ii + 1)*s1v(2, ik)
4756 ds1v(2, iat, 1, ik) = r(ii + 2)*s1v(1, ik)
4757 ds1v(2, iat, 2, ik) = r(ii + 2)*s1v(2, ik)
4758 ds1v(3, iat, 1, ik) = r(ii + 3)*s1v(1, ik)
4759 ds1v(3, iat, 2, ik) = r(ii + 3)*s1v(2, ik)
4766 DO iat = 1, rmsd_atom
4773 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4775 ALLOCATE (fi(3, rmsd_atom))
4777 DO iat = 1, rmsd_atom
4778 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))
4779 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))
4780 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))
4781 CALL put_derivative(colvar, iat, fi(:, iat))
4788 DEALLOCATE (vec_dif)
4789 DEALLOCATE (dvec_dif)
4794 END SUBROUTINE rpath_dist_rmsd
4801 SUBROUTINE rpath_rmsd(colvar, particles)
4805 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
4806 INTEGER,
DIMENSION(:),
POINTER :: iatom
4807 REAL(
dp) :: lambda, my_rmsd, s1(2)
4808 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0
4809 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: fi, riat, s1v
4810 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1
4811 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: ds1v
4812 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
4813 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
4814 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
4816 nconf = colvar%reaction_path_param%nr_frames
4817 rmsd_atom = colvar%reaction_path_param%n_components
4818 lambda = colvar%reaction_path_param%lambda
4819 path_conf => colvar%reaction_path_param%r_ref
4820 iatom => colvar%reaction_path_param%i_rmsd
4822 natom =
SIZE(particles)
4824 ALLOCATE (r0(3*natom))
4825 ALLOCATE (r(3*natom))
4826 ALLOCATE (riat(3, rmsd_atom))
4827 ALLOCATE (s1v(2, nconf))
4828 ALLOCATE (ds1v(3, rmsd_atom, 2, nconf))
4829 ALLOCATE (ds1(3, rmsd_atom, 2))
4830 ALLOCATE (drmsd(3, natom))
4832 ALLOCATE (weight(natom))
4836 r0(ii + 1) = particles(i)%r(1)
4837 r0(ii + 2) = particles(i)%r(2)
4838 r0(ii + 3) = particles(i)%r(3)
4841 DO iat = 1, rmsd_atom
4843 riat(:, iat) = particles(ii)%r
4848 DO iat = 1, rmsd_atom
4856 r(ii + 1) = path_conf(ii + 1, ik)
4857 r(ii + 2) = path_conf(ii + 2, ik)
4858 r(ii + 3) = path_conf(ii + 3, ik)
4861 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
4862 rotate=.false., drmsd3=drmsd)
4864 s1v(1, ik) = real(ik - 1,
dp)*exp(-lambda*my_rmsd)
4865 s1v(2, ik) = exp(-lambda*my_rmsd)
4866 DO iat = 1, rmsd_atom
4868 ds1v(1, iat, 1, ik) = drmsd(1, i)*s1v(1, ik)
4869 ds1v(1, iat, 2, ik) = drmsd(1, i)*s1v(2, ik)
4870 ds1v(2, iat, 1, ik) = drmsd(2, i)*s1v(1, ik)
4871 ds1v(2, iat, 2, ik) = drmsd(2, i)*s1v(2, ik)
4872 ds1v(3, iat, 1, ik) = drmsd(3, i)*s1v(1, ik)
4873 ds1v(3, iat, 2, ik) = drmsd(3, i)*s1v(2, ik)
4880 DO iat = 1, rmsd_atom
4887 colvar%ss = s1(1)/s1(2)/real(nconf - 1,
dp)
4889 ALLOCATE (fi(3, rmsd_atom))
4891 DO iat = 1, rmsd_atom
4892 fi(1, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(1, iat, 1) - ds1(1, iat, 2)*s1(1)/s1(2))
4893 fi(2, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(2, iat, 1) - ds1(2, iat, 2)*s1(1)/s1(2))
4894 fi(3, iat) = -lambda/s1(2)/real(nconf - 1,
dp)*(ds1(3, iat, 1) - ds1(3, iat, 2)*s1(1)/s1(2))
4895 CALL put_derivative(colvar, iat, fi(:, iat))
4908 END SUBROUTINE rpath_rmsd
4920 SUBROUTINE distance_from_path_colvar(colvar, cell, subsys, particles)
4925 POINTER :: particles
4931 IF (
PRESENT(particles))
THEN
4932 my_particles => particles
4934 cpassert(
PRESENT(subsys))
4936 my_particles => particles_i%els
4939 IF (colvar%reaction_path_param%dist_rmsd)
THEN
4940 CALL dpath_dist_rmsd(colvar, my_particles)
4941 ELSE IF (colvar%reaction_path_param%rmsd)
THEN
4942 CALL dpath_rmsd(colvar, my_particles)
4944 CALL dpath_colvar(colvar, cell, my_particles)
4947 END SUBROUTINE distance_from_path_colvar
4959 SUBROUTINE dpath_colvar(colvar, cell, particles)
4964 INTEGER :: i, iend, ii, istart, j, k, ncolv
4965 REAL(
dp) :: lambda, s1
4966 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: ds1, s1v, ss_vals
4967 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1v, f_vals, fi
4969 istart = colvar%reaction_path_param%function_bounds(1)
4970 iend = colvar%reaction_path_param%function_bounds(2)
4972 ncolv = colvar%reaction_path_param%n_components
4973 lambda = colvar%reaction_path_param%lambda
4974 ALLOCATE (f_vals(ncolv, istart:iend))
4975 f_vals(:, :) = colvar%reaction_path_param%f_vals
4976 ALLOCATE (ss_vals(ncolv))
4979 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar, cell, particles)
4980 ss_vals(i) = colvar%reaction_path_param%colvar_p(i)%colvar%ss
4983 ALLOCATE (s1v(istart:iend))
4984 ALLOCATE (ds1v(ncolv, istart:iend))
4985 ALLOCATE (ds1(ncolv))
4988 s1v(k) = exp(-lambda*dot_product(ss_vals(:) - f_vals(:, k), ss_vals(:) - f_vals(:, k)))
4990 ds1v(j, k) = f_vals(j, k)*s1v(k)
4998 colvar%ss = -1.0_dp/lambda*log(s1)
5000 ALLOCATE (fi(3, colvar%n_atom_s))
5004 DO j = 1, colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s
5006 fi(:, ii) = colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:, j)* &
5007 2.0_dp*(ss_vals(i) - ds1(i)/s1)
5011 DO i = 1, colvar%n_atom_s
5012 CALL put_derivative(colvar, i, fi(:, i))
5017 DEALLOCATE (ss_vals)
5022 END SUBROUTINE dpath_colvar
5033 SUBROUTINE dpath_dist_rmsd(colvar, particles)
5038 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5039 INTEGER,
DIMENSION(:),
POINTER :: iatom
5040 REAL(
dp) :: lambda, s1, sum_exp
5041 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v, vec_dif
5042 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, dvec_dif, fi, riat
5043 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5044 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5046 nconf = colvar%reaction_path_param%nr_frames
5047 rmsd_atom = colvar%reaction_path_param%n_components
5048 lambda = colvar%reaction_path_param%lambda
5049 path_conf => colvar%reaction_path_param%r_ref
5050 iatom => colvar%reaction_path_param%i_rmsd
5052 natom =
SIZE(particles)
5054 ALLOCATE (r0(3*natom))
5055 ALLOCATE (r(3*natom))
5056 ALLOCATE (riat(3, rmsd_atom))
5057 ALLOCATE (vec_dif(rmsd_atom))
5058 ALLOCATE (dvec_dif(3, rmsd_atom))
5059 ALLOCATE (s1v(nconf))
5060 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5061 ALLOCATE (ds1(3, rmsd_atom))
5064 r0(ii + 1) = particles(i)%r(1)
5065 r0(ii + 2) = particles(i)%r(2)
5066 r0(ii + 3) = particles(i)%r(3)
5069 DO iat = 1, rmsd_atom
5071 riat(:, iat) = particles(ii)%r
5077 r(ii + 1) = path_conf(ii + 1, ik)
5078 r(ii + 2) = path_conf(ii + 2, ik)
5079 r(ii + 3) = path_conf(ii + 3, ik)
5082 CALL rmsd3(particles, r, r0, output_unit=-1, rotate=.true.)
5085 DO iat = 1, rmsd_atom
5088 vec_dif(iat) = (riat(1, iat) - r(ii + 1))**2 + (riat(2, iat) - r(ii + 2))**2 + (riat(3, iat) - r(ii + 3))**2
5089 sum_exp = sum_exp + vec_dif(iat)
5090 dvec_dif(1, iat) = r(ii + 1)
5091 dvec_dif(2, iat) = r(ii + 2)
5092 dvec_dif(3, iat) = r(ii + 3)
5094 s1v(ik) = exp(-lambda*sum_exp)
5095 DO iat = 1, rmsd_atom
5096 ds1v(1, iat, ik) = dvec_dif(1, iat)*s1v(ik)
5097 ds1v(2, iat, ik) = dvec_dif(2, iat)*s1v(ik)
5098 ds1v(3, iat, ik) = dvec_dif(3, iat)*s1v(ik)
5103 DO iat = 1, rmsd_atom
5108 colvar%ss = -1.0_dp/lambda*log(s1)
5110 ALLOCATE (fi(3, rmsd_atom))
5112 DO iat = 1, rmsd_atom
5113 fi(:, iat) = 2.0_dp*(riat(:, iat) - ds1(:, iat)/s1)
5114 CALL put_derivative(colvar, iat, fi(:, iat))
5121 DEALLOCATE (vec_dif)
5122 DEALLOCATE (dvec_dif)
5126 END SUBROUTINE dpath_dist_rmsd
5133 SUBROUTINE dpath_rmsd(colvar, particles)
5138 INTEGER :: i, iat, ii, ik, natom, nconf, rmsd_atom
5139 INTEGER,
DIMENSION(:),
POINTER :: iatom
5140 REAL(
dp) :: lambda, my_rmsd, s1
5141 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: r, r0, s1v
5142 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ds1, fi, riat
5143 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ds1v
5144 REAL(
dp),
DIMENSION(:, :),
POINTER :: path_conf
5145 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weight
5146 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: drmsd
5148 nconf = colvar%reaction_path_param%nr_frames
5149 rmsd_atom = colvar%reaction_path_param%n_components
5150 lambda = colvar%reaction_path_param%lambda
5151 path_conf => colvar%reaction_path_param%r_ref
5152 iatom => colvar%reaction_path_param%i_rmsd
5154 natom =
SIZE(particles)
5156 ALLOCATE (r0(3*natom))
5157 ALLOCATE (r(3*natom))
5158 ALLOCATE (riat(3, rmsd_atom))
5159 ALLOCATE (s1v(nconf))
5160 ALLOCATE (ds1v(3, rmsd_atom, nconf))
5161 ALLOCATE (ds1(3, rmsd_atom))
5162 ALLOCATE (drmsd(3, natom))
5164 ALLOCATE (weight(natom))
5168 r0(ii + 1) = particles(i)%r(1)
5169 r0(ii + 2) = particles(i)%r(2)
5170 r0(ii + 3) = particles(i)%r(3)
5173 DO iat = 1, rmsd_atom
5175 riat(:, iat) = particles(ii)%r
5180 DO iat = 1, rmsd_atom
5188 r(ii + 1) = path_conf(ii + 1, ik)
5189 r(ii + 2) = path_conf(ii + 2, ik)
5190 r(ii + 3) = path_conf(ii + 3, ik)
5193 CALL rmsd3(particles, r0, r, output_unit=-1, weights=weight, my_val=my_rmsd, &
5194 rotate=.false., drmsd3=drmsd)
5196 s1v(ik) = exp(-lambda*my_rmsd)
5197 DO iat = 1, rmsd_atom
5199 ds1v(1, iat, ik) = drmsd(1, i)*s1v(ik)
5200 ds1v(2, iat, ik) = drmsd(2, i)*s1v(ik)
5201 ds1v(3, iat, ik) = drmsd(3, i)*s1v(ik)
5206 DO iat = 1, rmsd_atom
5211 colvar%ss = -1.0_dp/lambda*log(s1)
5213 ALLOCATE (fi(3, rmsd_atom))
5215 DO iat = 1, rmsd_atom
5216 fi(:, iat) = ds1(:, iat)/s1
5217 CALL put_derivative(colvar, iat, fi(:, iat))
5230 END SUBROUTINE dpath_rmsd
5241 SUBROUTINE population_colvar(colvar, cell, subsys, particles)
5246 POINTER :: particles
5248 INTEGER :: i, ii, jj, n_atoms_from, n_atoms_to, &
5250 REAL(
dp) :: dfunc, dfunc_coord, ftmp(3), func, func_coord, inv_n_atoms_from, invden, n_0, &
5251 ncoord, norm, num, population, r12, r_0, rdist, sigma, ss(3), xij(3)
5252 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftmp_coord
5253 REAL(
dp),
DIMENSION(3) :: xpi, xpj
5260 NULLIFY (particles_i)
5262 IF (
PRESENT(particles))
THEN
5263 my_particles => particles
5265 cpassert(
PRESENT(subsys))
5267 my_particles => particles_i%els
5269 n_atoms_to = colvar%population_param%n_atoms_to
5270 n_atoms_from = colvar%population_param%n_atoms_from
5271 nncrd = colvar%population_param%nncrd
5272 ndcrd = colvar%population_param%ndcrd
5273 r_0 = colvar%population_param%r_0
5274 n_0 = colvar%population_param%n0
5275 sigma = colvar%population_param%sigma
5277 ALLOCATE (ftmp_coord(3, n_atoms_to))
5283 colvar%dsdr = 0.0_dp
5284 inv_n_atoms_from = 1.0_dp/real(n_atoms_from, kind=
dp)
5286 norm = sqrt(
pi*2.0_dp)*sigma
5289 DO ii = 1, n_atoms_from
5290 i = colvar%population_param%i_at_from(ii)
5291 CALL get_coordinates(colvar, i, xpi, my_particles)
5292 DO jj = 1, n_atoms_to
5293 i = colvar%population_param%i_at_to(jj)
5294 CALL get_coordinates(colvar, i, xpj, my_particles)
5295 ss = matmul(cell%h_inv, xpi(:) - xpj(:))
5297 xij = matmul(cell%hmat, ss)
5298 r12 = sqrt(xij(1)**2 + xij(2)**2 + xij(3)**2)
5299 IF (r12 < 1.0e-8_dp) cycle
5301 num = (1.0_dp - rdist**nncrd)
5302 invden = 1.0_dp/(1.0_dp - rdist**ndcrd)
5303 func_coord = num*invden
5304 dfunc_coord = (-nncrd*rdist**(nncrd - 1)*invden &
5305 + num*(invden)**2*ndcrd*rdist**(ndcrd - 1))/(r12*r_0)
5307 ncoord = ncoord + func_coord
5308 ftmp_coord(1, jj) = dfunc_coord*xij(1)
5309 ftmp_coord(2, jj) = dfunc_coord*xij(2)
5310 ftmp_coord(3, jj) = dfunc_coord*xij(3)
5313 func = exp(-(ncoord - n_0)**2/(2.0_dp*sigma*sigma))
5314 dfunc = -func*(ncoord - n_0)/(sigma*sigma)
5316 population = population + norm*func
5317 DO jj = 1, n_atoms_to
5318 ftmp(1) = ftmp_coord(1, jj)*dfunc
5319 ftmp(2) = ftmp_coord(2, jj)*dfunc
5320 ftmp(3) = ftmp_coord(3, jj)*dfunc
5321 CALL put_derivative(colvar, ii, ftmp)
5322 ftmp(1) = -ftmp_coord(1, jj)*dfunc
5323 ftmp(2) = -ftmp_coord(2, jj)*dfunc
5324 ftmp(3) = -ftmp_coord(3, jj)*dfunc
5325 CALL put_derivative(colvar, n_atoms_from + jj, ftmp)
5329 colvar%ss = population
5330 END SUBROUTINE population_colvar
5342 SUBROUTINE gyration_radius_colvar(colvar, cell, subsys, particles)
5348 POINTER :: particles
5350 INTEGER :: i, ii, n_atoms
5351 REAL(
dp) :: dri2, func, gyration, inv_n, mass_tot, mi
5352 REAL(
dp),
DIMENSION(3) :: dfunc, dxi, ftmp, ss, xpcom, xpi
5356 NULLIFY (particles_i, my_particles)
5358 IF (
PRESENT(particles))
THEN
5359 my_particles => particles
5361 cpassert(
PRESENT(subsys))
5363 my_particles => particles_i%els
5365 n_atoms = colvar%gyration_param%n_atoms
5366 inv_n = 1.0_dp/n_atoms
5372 i = colvar%gyration_param%i_at(ii)
5373 CALL get_coordinates(colvar, i, xpi, my_particles)
5374 CALL get_mass(colvar, i, mi, my_particles)
5375 xpcom(:) = xpcom(:) + xpi(:)*mi
5376 mass_tot = mass_tot + mi
5378 xpcom(:) = xpcom(:)/mass_tot
5384 i = colvar%gyration_param%i_at(ii)
5385 CALL get_coordinates(colvar, i, xpi, my_particles)
5386 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5388 dxi = matmul(cell%hmat, ss)
5389 dri2 = (dxi(1)**2 + dxi(2)**2 + dxi(3)**2)
5391 dfunc(:) = dfunc(:) + dxi(:)
5393 gyration = sqrt(inv_n*func)
5396 i = colvar%gyration_param%i_at(ii)
5397 CALL get_coordinates(colvar, i, xpi, my_particles)
5398 CALL get_mass(colvar, i, mi, my_particles)
5399 ss = matmul(cell%h_inv, xpi(:) - xpcom(:))
5401 dxi = matmul(cell%hmat, ss)
5402 ftmp(1) = dxi(1) - dfunc(1)*mi/mass_tot
5403 ftmp(2) = dxi(2) - dfunc(2)*mi/mass_tot
5404 ftmp(3) = dxi(3) - dfunc(3)*mi/mass_tot
5405 ftmp(:) = ftmp(:)*inv_n/gyration
5406 CALL put_derivative(colvar, ii, ftmp)
5408 colvar%ss = gyration
5410 END SUBROUTINE gyration_radius_colvar
5421 SUBROUTINE rmsd_colvar(colvar, subsys, particles)
5425 POINTER :: particles
5427 CALL rmsd_colvar_low(colvar, subsys, particles)
5428 END SUBROUTINE rmsd_colvar
5443 SUBROUTINE rmsd_colvar_low(colvar, subsys, particles)
5448 POINTER :: particles
5450 INTEGER :: i, ii, natom, nframes
5451 REAL(kind=
dp) :: cv_val, f1, ftmp(3)
5452 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: der, r,
rmsd
5453 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r0
5454 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: drmsd
5455 REAL(kind=
dp),
DIMENSION(:),
POINTER :: weights
5459 NULLIFY (my_particles, particles_i, weights)
5461 IF (
PRESENT(particles))
THEN
5462 my_particles => particles
5464 cpassert(
PRESENT(subsys))
5466 my_particles => particles_i%els
5469 natom =
SIZE(my_particles)
5470 nframes = colvar%rmsd_param%nr_frames
5471 ALLOCATE (drmsd(3, natom, nframes))
5474 ALLOCATE (r0(3*natom, nframes))
5475 ALLOCATE (
rmsd(nframes))
5476 ALLOCATE (der(nframes))
5477 ALLOCATE (r(3*natom))
5479 weights => colvar%rmsd_param%weights
5482 r(ii + 1) = my_particles(i)%r(1)
5483 r(ii + 2) = my_particles(i)%r(2)
5484 r(ii + 3) = my_particles(i)%r(3)
5486 r0(:, :) = colvar%rmsd_param%r_ref
5489 CALL rmsd3(my_particles, r, r0(:, 1), output_unit=-1, weights=weights, my_val=
rmsd(1), rotate=.false., drmsd3=drmsd(:, :, 1))
5491 IF (nframes == 2)
THEN
5492 CALL rmsd3(my_particles, r, r0(:, 2), output_unit=-1, weights=weights, &
5493 my_val=
rmsd(2), rotate=.false., drmsd3=drmsd(:, :, 2))
5499 der(1) = f1 - cv_val*f1
5501 der(2) = -f1 - cv_val*f1
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) = der(1)*drmsd(1, ii, 1) + der(2)*drmsd(1, ii, 2)
5507 ftmp(2) = der(1)*drmsd(2, ii, 1) + der(2)*drmsd(2, ii, 2)
5508 ftmp(3) = der(1)*drmsd(3, ii, 1) + der(2)*drmsd(3, ii, 2)
5509 CALL put_derivative(colvar, i, ftmp)
5512 ELSE IF (nframes == 1)
THEN
5515 cv_val = sqrt(
rmsd(1))
5517 IF (cv_val /= 0.0_dp) f1 = 0.5_dp/cv_val
5518 DO i = 1, colvar%rmsd_param%n_atoms
5519 ii = colvar%rmsd_param%i_rmsd(i)
5520 IF (weights(ii) > 0.0_dp)
THEN
5521 ftmp(1) = f1*drmsd(1, ii, 1)
5522 ftmp(2) = f1*drmsd(2, ii, 1)
5523 ftmp(3) = f1*drmsd(3, ii, 1)
5524 CALL put_derivative(colvar, i, ftmp)
5528 cpabort(
"RMSD implemented only for 1 and 2 reference frames!")
5538 END SUBROUTINE rmsd_colvar_low
5550 SUBROUTINE ring_puckering_colvar(colvar, cell, subsys, particles)
5555 POINTER :: particles
5557 INTEGER :: i, ii, j, jj, m, nring
5558 REAL(kind=
dp) :: a, at, b, da, db, ds, kr, rpxpp, svar
5559 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: cosj, sinj, z
5560 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: r
5561 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: nforce, zforce
5562 REAL(kind=
dp),
DIMENSION(3) :: ftmp, nv, r0, rp, rpp, uv
5563 REAL(kind=
dp),
DIMENSION(3, 3) :: dnvp, dnvpp
5568 IF (
PRESENT(particles))
THEN
5569 my_particles => particles
5571 cpassert(
PRESENT(subsys))
5573 my_particles => particles_i%els
5576 nring = colvar%ring_puckering_param%nring
5577 ALLOCATE (r(3, nring), z(nring), cosj(nring), sinj(nring))
5578 ALLOCATE (nforce(3, 3, nring), zforce(nring, nring, 3))
5580 i = colvar%ring_puckering_param%atoms(ii)
5581 CALL get_coordinates(colvar, i, r(:, ii), my_particles)
5586 r(:, ii) =
pbc(r(:, ii), r0, cell)
5591 r0(:) = r0(:) + r(:, ii)
5593 kr = 1._dp/real(nring, kind=
dp)
5596 r(:, ii) = r(:, ii) - r0(:)
5602 cosj(ii) = cos(
twopi*(ii - 1)*kr)
5603 sinj(ii) = sin(
twopi*(ii - 1)*kr)
5604 rp(:) = rp(:) + r(:, ii)*sinj(ii)
5605 rpp(:) = rpp(:) + r(:, ii)*cosj(ii)
5608 nv = nv/sqrt(sum(nv**2))
5612 rpxpp = sqrt(sum(uv**2))
5617 dnvp(:, i) = uv - nv*sum(uv*nv)
5621 dnvpp(:, i) = uv - nv*sum(uv*nv)
5624 nforce(:, :, ii) = dnvp(:, :)*sinj(ii) + dnvpp(:, :)*cosj(ii)
5629 z(ii) = sum(r(:, ii)*nv(:))
5635 zforce(ii, jj, :) = nv
5637 zforce(ii, jj, :) = 0._dp
5641 zforce(ii, jj, i) = zforce(ii, jj, i) + r(j, ii)*nforce(j, i, jj)
5647 IF (colvar%ring_puckering_param%iq == 0)
THEN
5649 svar = sqrt(sum(z**2))
5653 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*z(jj)
5656 CALL put_derivative(colvar, ii, ftmp)
5659 m = abs(colvar%ring_puckering_param%iq)
5661 IF (mod(nring, 2) == 0 .AND. colvar%ring_puckering_param%iq == nring/2)
THEN
5665 IF (mod(ii, 2) == 0)
THEN
5671 svar = svar*sqrt(kr)
5675 IF (mod(jj, 2) == 0)
THEN
5676 ftmp(:) = ftmp(:) - zforce(jj, ii, :)*sqrt(kr)
5678 ftmp(:) = ftmp(:) + zforce(jj, ii, :)*sqrt(kr)
5681 CALL put_derivative(colvar, ii, -ftmp)
5684 cpassert(m <= (nring - 1)/2)
5688 a = a + z(ii)*cos(
twopi*m*(ii - 1)*kr)
5689 b = b - z(ii)*sin(
twopi*m*(ii - 1)*kr)
5691 a = a*sqrt(2._dp*kr)
5692 b = b*sqrt(2._dp*kr)
5693 IF (colvar%ring_puckering_param%iq > 0)
THEN
5695 svar = sqrt(a*a + b*b)
5701 IF (at >
pi/2._dp)
THEN
5702 svar = 2.5_dp*
pi - at
5704 svar = 0.5_dp*
pi - at
5712 ds = da*cos(
twopi*m*(ii - 1)*kr)
5713 ds = ds - db*sin(
twopi*m*(ii - 1)*kr)
5714 ftmp(:) = ftmp(:) + ds*sqrt(2._dp*kr)*zforce(ii, jj, :)
5716 CALL put_derivative(colvar, jj, ftmp)
5723 DEALLOCATE (r, z, cosj, sinj, nforce, zforce)
5725 END SUBROUTINE ring_puckering_colvar
5747 RECURSIVE FUNCTION rec_eval_grid(iw1, ncol, f_vals, v_count, &
5748 gp, grid_sp, step_size, istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
RESULT(k)
5749 INTEGER :: iw1, ncol
5750 REAL(
dp),
DIMENSION(:, :),
POINTER :: f_vals
5752 REAL(
dp),
DIMENSION(:),
POINTER :: gp, grid_sp
5753 REAL(
dp) :: step_size
5754 INTEGER :: istart, iend
5755 REAL(
dp),
DIMENSION(:, :),
POINTER :: s1v
5756 REAL(
dp),
DIMENSION(:),
POINTER :: s1
5757 INTEGER,
DIMENSION(:, :),
POINTER :: p_bounds
5759 INTEGER :: ifunc, nconf, k
5761 INTEGER :: count1, i
5764 IF (v_count < ncol)
THEN
5765 count1 = v_count + 1
5766 DO i = p_bounds(1, count1), p_bounds(2, count1)
5767 gp(count1) = real(i, kind=
dp)*grid_sp(count1)
5768 k = rec_eval_grid(iw1, ncol, f_vals, count1, gp, grid_sp, step_size, &
5769 istart, iend, s1v, s1, p_bounds, lambda, ifunc, nconf)
5771 ELSE IF (v_count == ncol .AND. ifunc == 1)
THEN
5773 s1v(1, i) = real(i, kind=
dp)*step_size*exp(-lambda*dot_product(gp(:) - f_vals(:, i), &
5774 gp(:) - f_vals(:, i)))
5775 s1v(2, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5780 WRITE (iw1,
'(5F10.5)') gp(:), s1(1)/s1(2)/real(nconf - 1,
dp)
5781 ELSE IF (v_count == ncol .AND. ifunc == 2)
THEN
5783 s1v(1, i) = exp(-lambda*dot_product(gp(:) - f_vals(:, i), gp(:) - f_vals(:, i)))
5787 WRITE (iw1,
'(5F10.5)') gp(:), -lambda*log(s1(1))
5789 END FUNCTION rec_eval_grid
5802 SUBROUTINE read_frames(frame_section, para_env, nr_frames, r_ref, n_atoms)
5806 INTEGER,
INTENT(IN) :: nr_frames
5807 REAL(
dp),
DIMENSION(:, :),
POINTER :: r_ref
5808 INTEGER,
INTENT(OUT) :: n_atoms
5810 CHARACTER(LEN=default_path_length) :: filename
5811 CHARACTER(LEN=default_string_length) :: dummy_char
5812 INTEGER :: i, j, natom
5813 LOGICAL :: explicit, my_end
5814 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rptr
5827 ALLOCATE (r_ref(3*natom, nr_frames))
5830 cpassert(3*natom ==
SIZE(r_ref, 1))
5834 i_rep_val=j, r_vals=rptr)
5835 r_ref((j - 1)*3 + 1:(j - 1)*3 + 3, i) = rptr(1:3)
5841 cpassert(trim(filename) /=
"")
5843 CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
5849 ALLOCATE (r_ref(3*natom, nr_frames))
5852 cpassert(3*natom ==
SIZE(r_ref, 1))
5858 CALL cp_abort(__location__, &
5859 "Number of lines in XYZ format not equal to the number of atoms."// &
5860 " Error in XYZ format for COORD_A (CV rmsd). Very probably the"// &
5861 " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
5863 READ (parser%input_line, *) dummy_char, rptr(1:3)
5874 END SUBROUTINE read_frames
5885 SUBROUTINE wc_colvar(colvar, cell, subsys, particles, qs_env)
5886 TYPE(colvar_type),
POINTER :: colvar
5887 TYPE(cell_type),
POINTER :: cell
5888 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5889 TYPE(particle_type),
DIMENSION(:), &
5890 OPTIONAL,
POINTER :: particles
5891 TYPE(qs_environment_type),
POINTER,
OPTIONAL :: qs_env
5893 INTEGER :: od, h, oa
5894 REAL(dp) :: rod(3), roa(3), rh(3), &
5895 x, y, s(3), xv(3), dmin, amin
5896 INTEGER :: idmin, iamin, i, j
5897 TYPE(particle_list_type),
POINTER :: particles_i
5898 TYPE(particle_type),
DIMENSION(:), &
5899 POINTER :: my_particles
5900 TYPE(wannier_centres_type),
DIMENSION(:),
POINTER :: wc
5901 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
5902 INTEGER :: nwca, nwcd
5905 NULLIFY (particles_i, wc)
5907 cpassert(colvar%type_id == wc_colvar_id)
5908 IF (
PRESENT(particles))
THEN
5909 my_particles => particles
5911 cpassert(
PRESENT(subsys))
5912 CALL cp_subsys_get(subsys, particles=particles_i)
5913 my_particles => particles_i%els
5915 CALL get_qs_env(qs_env, wanniercentres=wc)
5916 rcut = colvar%Wc%rcut
5917 od = colvar%Wc%ids(1)
5918 h = colvar%Wc%ids(2)
5919 oa = colvar%Wc%ids(3)
5920 CALL get_coordinates(colvar, od, rod, my_particles)
5921 CALL get_coordinates(colvar, h, rh, my_particles)
5922 CALL get_coordinates(colvar, oa, roa, my_particles)
5923 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
5924 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
5927 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
5928 x = distance(rod - wc(1)%centres(:, j))
5929 y = distance(roa - wc(1)%centres(:, j))
5941 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
5942 amin = distance(rh - wc(1)%centres(:, wcai(1)))
5947 x = distance(rh - wc(1)%centres(:, wcdi(i)))
5954 x = distance(rh - wc(1)%centres(:, wcai(i)))
5961 colvar%ss = wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
5971 REAL(dp) function distance(rij)
5972 REAL(dp),
INTENT(in) :: rij(3)
5974 s = matmul(cell%h_inv, rij)
5976 xv = matmul(cell%hmat, s)
5977 distance = norm2(xv)
5978 END FUNCTION distance
5980 END SUBROUTINE wc_colvar
5991 SUBROUTINE hbp_colvar(colvar, cell, subsys, particles, qs_env)
5992 TYPE(colvar_type),
POINTER :: colvar
5993 TYPE(cell_type),
POINTER :: cell
5994 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
5995 TYPE(particle_type),
DIMENSION(:), &
5996 OPTIONAL,
POINTER :: particles
5997 TYPE(qs_environment_type),
OPTIONAL,
POINTER :: qs_env
5999 INTEGER :: od, h, oa
6000 REAL(dp) :: rod(3), roa(3), rh(3), &
6001 x, y, s(3), xv(3), dmin, amin
6002 INTEGER :: idmin, iamin, i, j, il, output_unit
6003 TYPE(particle_list_type),
POINTER :: particles_i
6004 TYPE(particle_type),
DIMENSION(:), &
6005 POINTER :: my_particles
6006 TYPE(wannier_centres_type), &
6007 DIMENSION(:),
POINTER :: wc
6008 INTEGER,
ALLOCATABLE :: wcai(:), wcdi(:)
6009 INTEGER :: nwca, nwcd
6012 NULLIFY (particles_i, wc)
6013 output_unit = cp_logger_get_default_io_unit()
6015 cpassert(colvar%type_id == hbp_colvar_id)
6016 IF (
PRESENT(particles))
THEN
6017 my_particles => particles
6019 cpassert(
PRESENT(subsys))
6020 CALL cp_subsys_get(subsys, particles=particles_i)
6021 my_particles => particles_i%els
6023 CALL get_qs_env(qs_env, wanniercentres=wc)
6024 rcut = colvar%HBP%rcut
6025 ALLOCATE (wcai(
SIZE(wc(1)%WannierHamDiag)))
6026 ALLOCATE (wcdi(
SIZE(wc(1)%WannierHamDiag)))
6028 DO il = 1, colvar%HBP%nPoints
6029 od = colvar%HBP%ids(il, 1)
6030 h = colvar%HBP%ids(il, 2)
6031 oa = colvar%HBP%ids(il, 3)
6032 CALL get_coordinates(colvar, od, rod, my_particles)
6033 CALL get_coordinates(colvar, h, rh, my_particles)
6034 CALL get_coordinates(colvar, oa, roa, my_particles)
6037 DO j = 1,
SIZE(wc(1)%WannierHamDiag)
6038 x = distance(rod - wc(1)%centres(:, j))
6039 y = distance(roa - wc(1)%centres(:, j))
6051 dmin = distance(rh - wc(1)%centres(:, wcdi(1)))
6052 amin = distance(rh - wc(1)%centres(:, wcai(1)))
6057 x = distance(rh - wc(1)%centres(:, wcdi(i)))
6064 x = distance(rh - wc(1)%centres(:, wcai(i)))
6070 colvar%HBP%ewc(il) = colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6071 colvar%ss = colvar%ss + colvar%HBP%shift + wc(1)%WannierHamDiag(idmin) - wc(1)%WannierHamDiag(iamin)
6073 IF (output_unit > 0)
THEN
6074 DO il = 1, colvar%HBP%nPoints
6075 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP| = ", colvar%HBP%ewc(il)
6077 WRITE (output_unit,
'(a,1(f16.8,1x))')
"HBP|\theta(x) = ", colvar%ss
6088 REAL(dp) function distance(rij)
6089 REAL(dp),
INTENT(in) :: rij(3)
6091 s = matmul(cell%h_inv, rij)
6093 xv = matmul(cell%hmat, s)
6094 distance = norm2(xv)
6095 END FUNCTION distance
6097 END SUBROUTINE hbp_colvar
6106 SUBROUTINE colvar_eval_voronoiipz(colvar, cell, subsys, particles)
6107 TYPE(colvar_type),
POINTER :: colvar
6108 TYPE(cell_type),
POINTER :: cell
6109 TYPE(cp_subsys_type),
OPTIONAL,
POINTER :: subsys
6110 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
6111 POINTER :: particles
6113 INTEGER :: i, iatom, j, jatom, n_atoms_a, &
6114 n_atoms_o, n_atoms_b, zidx
6115 REAL(kind=dp) :: charge_deriv, cutoff, denominator, &
6116 distance, exponent_max, factor, ion_z, &
6117 lambda, weighted_charge_deriv, zdist, zmid
6118 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: charge, charge_factor
6119 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: inv_distance, weight
6120 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: distance_vector
6121 TYPE(particle_type),
DIMENSION(:),
POINTER :: my_particles
6123 cpassert(colvar%type_id == voronoiipz_colvar_id)
6124 IF (
PRESENT(particles))
THEN
6125 my_particles => particles
6127 cpassert(
PRESENT(subsys))
6128 my_particles => subsys%particles%els
6131 n_atoms_a =
SIZE(colvar%voronoiipz_params%group_a)
6132 n_atoms_b =
SIZE(colvar%voronoiipz_params%group_b)
6133 n_atoms_o = n_atoms_a - colvar%voronoiipz_params%nrx
6134 cutoff = colvar%voronoiipz_params%nl_cutoff
6135 lambda = colvar%voronoiipz_params%lambda
6136 zidx = colvar%voronoiipz_params%zidx
6137 zmid = colvar%voronoiipz_params%zmid
6139 ALLOCATE (charge(n_atoms_a), charge_factor(n_atoms_a))
6140 ALLOCATE (distance_vector(3, n_atoms_a, n_atoms_b))
6141 ALLOCATE (inv_distance(n_atoms_a, n_atoms_b), weight(n_atoms_a, n_atoms_b))
6143 charge_factor = 0.0_dp
6144 distance_vector = 0.0_dp
6145 inv_distance = 0.0_dp
6151 jatom = colvar%voronoiipz_params%group_b(j)
6152 exponent_max = -huge(0.0_dp)
6154 iatom = colvar%voronoiipz_params%group_a(i)
6155 distance_vector(:, i, j) = pbc(my_particles(iatom)%r, my_particles(jatom)%r, cell)
6156 distance = norm2(distance_vector(:, i, j))
6157 IF (distance <= cutoff)
THEN
6158 IF (distance <= epsilon(0.0_dp))
THEN
6159 cpabort(
"VORONOIIPZ is undefined for coincident GROUPA and GROUPB atoms")
6161 inv_distance(i, j) = 1.0_dp/distance
6162 weight(i, j) = lambda*distance
6163 exponent_max = max(exponent_max, weight(i, j))
6167 denominator = 0.0_dp
6169 IF (inv_distance(i, j) > 0.0_dp)
THEN
6170 weight(i, j) = exp(weight(i, j) - exponent_max)
6171 denominator = denominator + weight(i, j)
6174 IF (denominator <= 0.0_dp)
THEN
6175 cpabort(
"VORONOIIPZ found a GROUPB atom without a GROUPA neighbor inside NL_CUTOFF")
6177 weight(:, j) = weight(:, j)/denominator
6180 charge(:) = sum(weight, dim=2)
6182 IF (i <= n_atoms_o)
THEN
6183 charge(i) = charge(i) - colvar%voronoiipz_params%d0
6185 SELECT CASE (i - n_atoms_o)
6187 charge(i) = charge(i) - colvar%voronoiipz_params%d1
6189 charge(i) = charge(i) - colvar%voronoiipz_params%d2
6191 charge(i) = charge(i) - colvar%voronoiipz_params%d3
6198 IF (charge(i) > 0.0_dp)
THEN
6199 iatom = colvar%voronoiipz_params%group_a(i)
6200 zdist = my_particles(iatom)%r(zidx) - zmid
6201 ion_z = ion_z + abs(zdist)*charge(i)**2
6202 charge_factor(i) = 2.0_dp*abs(zdist)*charge(i)
6203 IF (zdist > 0.0_dp)
THEN
6204 colvar%dsdr(zidx, i) = colvar%dsdr(zidx, i) + charge(i)**2
6205 ELSE IF (zdist < 0.0_dp)
THEN
6206 colvar%dsdr(zidx, i) = colvar%dsdr(zidx, i) - charge(i)**2
6213 weighted_charge_deriv = sum(weight(:, j)*charge_factor)
6215 IF (weight(i, j) > 0.0_dp)
THEN
6216 charge_deriv = charge_factor(i) - weighted_charge_deriv
6217 factor = lambda*weight(i, j)*charge_deriv*inv_distance(i, j)
6218 colvar%dsdr(:, i) = colvar%dsdr(:, i) - factor*distance_vector(:, i, j)
6219 colvar%dsdr(:, n_atoms_a + j) = colvar%dsdr(:, n_atoms_a + j) + &
6220 factor*distance_vector(:, i, j)
6226 DEALLOCATE (charge, charge_factor, distance_vector, inv_distance, weight)
6228 END SUBROUTINE colvar_eval_voronoiipz
6235 SUBROUTINE read_voronoiipz_colvars(voronoiipz_section, colvar)
6236 TYPE(section_vals_type),
POINTER :: voronoiipz_section
6237 TYPE(colvar_type),
POINTER :: colvar
6239 INTEGER :: i, j, n_atoms_a
6240 INTEGER,
DIMENSION(:),
POINTER :: atom_list
6242 CALL section_vals_val_get(voronoiipz_section,
"LAMBDA", r_val=colvar%voronoiipz_params%lambda)
6243 CALL section_vals_val_get(voronoiipz_section,
"ZIDX", i_val=colvar%voronoiipz_params%zidx)
6244 CALL section_vals_val_get(voronoiipz_section,
"NRX", i_val=colvar%voronoiipz_params%nrx)
6245 CALL section_vals_val_get(voronoiipz_section,
"ZMID", r_val=colvar%voronoiipz_params%zmid)
6246 CALL section_vals_val_get(voronoiipz_section,
"D_0", r_val=colvar%voronoiipz_params%d0)
6247 CALL section_vals_val_get(voronoiipz_section,
"D_1", r_val=colvar%voronoiipz_params%d1)
6248 CALL section_vals_val_get(voronoiipz_section,
"D_2", r_val=colvar%voronoiipz_params%d2)
6249 CALL section_vals_val_get(voronoiipz_section,
"D_3", r_val=colvar%voronoiipz_params%d3)
6250 CALL section_vals_val_get(voronoiipz_section,
"NL_CUTOFF", r_val=colvar%voronoiipz_params%nl_cutoff)
6253 CALL section_vals_val_get(voronoiipz_section,
"GROUPA", i_vals=atom_list)
6254 IF (.NOT.
ASSOCIATED(atom_list))
THEN
6255 cpabort(
"VORONOIIPZ requires a non-empty GROUPA")
6257 IF (
SIZE(atom_list) == 0) cpabort(
"VORONOIIPZ requires a non-empty GROUPA")
6258 ALLOCATE (colvar%voronoiipz_params%group_a(
SIZE(atom_list)))
6259 colvar%voronoiipz_params%group_a = atom_list
6262 CALL section_vals_val_get(voronoiipz_section,
"GROUPB", i_vals=atom_list)
6263 IF (.NOT.
ASSOCIATED(atom_list))
THEN
6264 cpabort(
"VORONOIIPZ requires a non-empty GROUPB")
6266 IF (
SIZE(atom_list) == 0) cpabort(
"VORONOIIPZ requires a non-empty GROUPB")
6267 ALLOCATE (colvar%voronoiipz_params%group_b(
SIZE(atom_list)))
6268 colvar%voronoiipz_params%group_b = atom_list
6270 n_atoms_a =
SIZE(colvar%voronoiipz_params%group_a)
6271 IF (colvar%voronoiipz_params%zidx < 1 .OR. colvar%voronoiipz_params%zidx > 3)
THEN
6272 cpabort(
"VORONOIIPZ ZIDX must be 1, 2, or 3")
6274 IF (colvar%voronoiipz_params%nrx < 0 .OR. colvar%voronoiipz_params%nrx > 3 .OR. &
6275 colvar%voronoiipz_params%nrx >= n_atoms_a)
THEN
6276 cpabort(
"VORONOIIPZ NRX must be between zero and three and smaller than the size of GROUPA")
6278 IF (colvar%voronoiipz_params%nl_cutoff <= 0.0_dp)
THEN
6279 cpabort(
"VORONOIIPZ NL_CUTOFF must be positive")
6281 IF (any(colvar%voronoiipz_params%group_a <= 0) .OR. &
6282 any(colvar%voronoiipz_params%group_b <= 0))
THEN
6283 cpabort(
"VORONOIIPZ atom indices must be positive")
6286 IF (any(colvar%voronoiipz_params%group_a(i) == colvar%voronoiipz_params%group_a(i + 1:)))
THEN
6287 cpabort(
"VORONOIIPZ GROUPA contains duplicate atom indices")
6289 IF (any(colvar%voronoiipz_params%group_a(i) == colvar%voronoiipz_params%group_b))
THEN
6290 cpabort(
"VORONOIIPZ GROUPA and GROUPB must be disjoint")
6293 DO j = 1,
SIZE(colvar%voronoiipz_params%group_b)
6294 IF (any(colvar%voronoiipz_params%group_b(j) == colvar%voronoiipz_params%group_b(j + 1:)))
THEN
6295 cpabort(
"VORONOIIPZ GROUPB contains duplicate atom indices")
6299 END SUBROUTINE read_voronoiipz_colvars
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 voronoiipz_colvar_id
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
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, matrix_vhxc, 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
subroutine check(condition, label)
Stops with an error message if a test fails.
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