|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00014 MODULE colvar_methods 00015 USE cell_types, ONLY: cell_type,& 00016 pbc 00017 USE colvar_types, ONLY: & 00018 HBP_colvar_id, Wc_colvar_id, angle_colvar_id, colvar_check_points, & 00019 colvar_create, colvar_setup, colvar_type, combine_colvar_id, & 00020 coord_colvar_id, dfunct_colvar_id, dist_colvar_id, & 00021 distance_from_path_colvar_id, eval_point_der, eval_point_mass, & 00022 eval_point_pos, gyration_colvar_id, hydronium_colvar_id, & 00023 mindist_colvar_id, plane_distance_colvar_id, & 00024 plane_plane_angle_colvar_id, population_colvar_id, qparm_colvar_id, & 00025 reaction_path_colvar_id, ring_puckering_colvar_id, rmsd_colvar_id, & 00026 rotation_colvar_id, torsion_colvar_id, u_colvar_id, & 00027 xyz_diag_colvar_id, xyz_outerdiag_colvar_id 00028 USE constraint_fxd, ONLY: check_fixed_atom_cns_colv 00029 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00030 cp_print_key_unit_nr 00031 USE cp_para_types, ONLY: cp_para_env_type 00032 USE cp_parser_methods, ONLY: parser_get_next_line,& 00033 parser_get_object 00034 USE cp_parser_types, ONLY: cp_parser_type,& 00035 parser_create,& 00036 parser_release 00037 USE cp_subsys_types, ONLY: cp_subsys_get,& 00038 cp_subsys_p_type,& 00039 cp_subsys_type 00040 USE cp_units, ONLY: cp_unit_to_cp2k 00041 USE f77_blas 00042 USE force_env_types, ONLY: force_env_get,& 00043 force_env_type,& 00044 use_mixed_force 00045 USE force_fields_util, ONLY: get_generic_info 00046 USE fparser, ONLY: EvalErrType,& 00047 evalf,& 00048 evalfd,& 00049 finalizef,& 00050 initf,& 00051 parsef 00052 USE input_constants, ONLY: & 00053 do_clv_x, do_clv_xy, do_clv_xz, do_clv_y, do_clv_yz, do_clv_z, & 00054 plane_def_atoms, plane_def_vec, rmsd_all, rmsd_list, rmsd_weightlist 00055 USE input_cp2k_colvar, ONLY: create_colvar_xyz_d_section,& 00056 create_colvar_xyz_od_section 00057 USE input_enumeration_types, ONLY: enum_i2c,& 00058 enumeration_type 00059 USE input_keyword_types, ONLY: keyword_get,& 00060 keyword_type 00061 USE input_section_types, ONLY: section_get_keyword,& 00062 section_release,& 00063 section_type,& 00064 section_vals_get,& 00065 section_vals_get_subs_vals,& 00066 section_vals_type,& 00067 section_vals_val_get 00068 USE kahan_sum, ONLY: accurate_sum 00069 USE kinds, ONLY: default_path_length,& 00070 default_string_length,& 00071 dp 00072 USE mathconstants, ONLY: fac,& 00073 maxfac,& 00074 pi,& 00075 twopi 00076 USE mathlib, ONLY: vector_product 00077 USE memory_utilities, ONLY: reallocate 00078 USE message_passing, ONLY: mp_sum,& 00079 mp_sync 00080 USE mixed_energy_types, ONLY: mixed_force_type 00081 USE mixed_environment_utils, ONLY: get_subsys_map_index 00082 USE molecule_kind_types, ONLY: fixd_constraint_type 00083 USE particle_list_types, ONLY: particle_list_p_type,& 00084 particle_list_type 00085 USE particle_types, ONLY: particle_type 00086 USE qs_environment_types, ONLY: get_qs_env,& 00087 qs_environment_type 00088 USE rmsd, ONLY: rmsd3 00089 USE spherical_harmonics, ONLY: dlegendre,& 00090 legendre 00091 USE string_utilities, ONLY: compress,& 00092 uppercase 00093 USE termination, ONLY: stop_program 00094 USE timings, ONLY: timeset,& 00095 timestop 00096 USE wannier_states_types, ONLY: wannier_centres_type 00097 #include "cp_common_uses.h" 00098 00099 IMPLICIT NONE 00100 PRIVATE 00101 00102 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'colvar_methods' 00103 REAL(KIND=dp), PRIVATE, PARAMETER :: tolerance_acos = 1.0E-5_dp 00104 00105 PUBLIC :: colvar_read,& 00106 colvar_eval_glob_f,& 00107 colvar_eval_mol_f 00108 00109 CONTAINS 00110 00111 ! ***************************************************************************** 00124 RECURSIVE SUBROUTINE colvar_read(colvar, icol, colvar_section, para_env, error) 00125 TYPE(colvar_type), POINTER :: colvar 00126 INTEGER, INTENT(IN) :: icol 00127 TYPE(section_vals_type), POINTER :: colvar_section 00128 TYPE(cp_para_env_type), POINTER :: para_env 00129 TYPE(cp_error_type), INTENT(inout) :: error 00130 00131 CHARACTER(len=*), PARAMETER :: routineN = 'colvar_read', 00132 routineP = moduleN//':'//routineN 00133 00134 CHARACTER(LEN=3) :: fmid 00135 CHARACTER(LEN=7) :: tag, tag_comp, tag_comp1, 00136 tag_comp2 00137 CHARACTER(LEN=default_path_length) :: path_function 00138 CHARACTER(LEN=default_string_length) :: tmpStr, tmpStr2 00139 CHARACTER(LEN=default_string_length), 00140 DIMENSION(:), POINTER :: c_kinds, my_par 00141 INTEGER :: handle, i, iatm, icomponent, iend, ifunc, ii, isize, istart, 00142 iw, iw1, j, k, kk, n_var, n_var_k, ncol, ndim, nr_frame, stat, v_count 00143 INTEGER, DIMENSION(:), POINTER :: iatms 00144 INTEGER, DIMENSION(:, :), POINTER :: p_bounds 00145 LOGICAL :: check, failure, 00146 use_mixed_energy 00147 LOGICAL, DIMENSION(23) :: my_subsection 00148 REAL(dp), DIMENSION(:), POINTER :: s1, wei, weights 00149 REAL(dp), DIMENSION(:, :), POINTER :: p_range, s1v 00150 REAL(KIND=dp), DIMENSION(1) :: my_val 00151 REAL(KIND=dp), DIMENSION(:), POINTER :: g_range, grid_point, grid_sp, 00152 my_vals, range 00153 TYPE(cp_logger_type), POINTER :: logger 00154 TYPE(enumeration_type), POINTER :: enum 00155 TYPE(keyword_type), POINTER :: keyword 00156 TYPE(section_type), POINTER :: section 00157 TYPE(section_vals_type), POINTER :: angle_section, colvar_subsection, 00158 combine_section, coordination_section, dfunct_section, 00159 distance_from_path_section, distance_section, frame_section, 00160 gyration_section, HBP_section, hydronium_section, mindist_section, 00161 path_section, plane_dist_section, plane_plane_angle_section, 00162 plane_sections, point_section, population_section, qparm_section, 00163 reaction_path_section, ring_puckering_section, rmsd_section, 00164 rotation_section, torsion_section, u_section, Wc_section, wrk_section, 00165 xyz_diag_section, xyz_outerdiag_section 00166 00167 failure=.FALSE. 00168 CALL timeset(routineN,handle) 00169 NULLIFY(logger, c_kinds, iatms) 00170 logger => cp_error_get_logger(error) 00171 my_subsection = .FALSE. 00172 distance_section => section_vals_get_subs_vals(colvar_section,"DISTANCE",i_rep_section=icol,error=error) 00173 dfunct_section => section_vals_get_subs_vals(colvar_section,"DISTANCE_FUNCTION",& 00174 i_rep_section=icol,error=error) 00175 angle_section => section_vals_get_subs_vals(colvar_section,"ANGLE",i_rep_section=icol,error=error) 00176 torsion_section => section_vals_get_subs_vals(colvar_section,"TORSION",i_rep_section=icol,error=error) 00177 coordination_section => section_vals_get_subs_vals(colvar_section,"COORDINATION",i_rep_section=icol,& 00178 error=error) 00179 plane_dist_section => section_vals_get_subs_vals(colvar_section,"DISTANCE_POINT_PLANE",i_rep_section=icol,& 00180 error=error) 00181 plane_plane_angle_section & 00182 => section_vals_get_subs_vals(colvar_section,"ANGLE_PLANE_PLANE",i_rep_section=icol,& 00183 error=error) 00184 rotation_section => section_vals_get_subs_vals(colvar_section,"BOND_ROTATION",i_rep_section=icol,& 00185 error=error) 00186 qparm_section => section_vals_get_subs_vals(colvar_section,"QPARM",i_rep_section=icol,error=error) 00187 hydronium_section => section_vals_get_subs_vals(colvar_section,"HYDRONIUM",i_rep_section=icol,error=error) 00188 reaction_path_section => section_vals_get_subs_vals(colvar_section,"REACTION_PATH",i_rep_section=icol,& 00189 can_return_null=.TRUE.,error=error) 00190 distance_from_path_section & 00191 => section_vals_get_subs_vals(colvar_section,"DISTANCE_FROM_PATH",& 00192 i_rep_section=icol, can_return_null=.TRUE.,error=error) 00193 combine_section => section_vals_get_subs_vals(colvar_section,"COMBINE_COLVAR",i_rep_section=icol,& 00194 can_return_null=.TRUE.,error=error) 00195 population_section => section_vals_get_subs_vals(colvar_section,"POPULATION",i_rep_section=icol,error=error) 00196 gyration_section => section_vals_get_subs_vals(colvar_section,"GYRATION_RADIUS",i_rep_section=icol,& 00197 error=error) 00198 rmsd_section => section_vals_get_subs_vals(colvar_section,"RMSD",i_rep_section=icol,error=error) 00199 xyz_diag_section => section_vals_get_subs_vals(colvar_section,"XYZ_DIAG",i_rep_section=icol,error=error) 00200 xyz_outerdiag_section => section_vals_get_subs_vals(colvar_section,"XYZ_OUTERDIAG",i_rep_section=icol,& 00201 error=error) 00202 u_section => section_vals_get_subs_vals(colvar_section,"U",i_rep_section=icol,error=error) 00203 Wc_section => section_vals_get_subs_vals(colvar_section,"WC",i_rep_section=icol,error=error) 00204 HBP_section => section_vals_get_subs_vals(colvar_section,"HBP",i_rep_section=icol,error=error) 00205 ring_puckering_section& 00206 => section_vals_get_subs_vals(colvar_section,"RING_PUCKERING",i_rep_section=icol,error=error) 00207 mindist_section => section_vals_get_subs_vals(colvar_section,"CONDITIONED_DISTANCE",i_rep_section=icol,error=error) 00208 00209 CALL section_vals_get(distance_section, explicit=my_subsection( 1), error=error) 00210 CALL section_vals_get(angle_section, explicit=my_subsection( 2), error=error) 00211 CALL section_vals_get(torsion_section, explicit=my_subsection( 3), error=error) 00212 CALL section_vals_get(coordination_section, explicit=my_subsection( 4), error=error) 00213 CALL section_vals_get(plane_dist_section, explicit=my_subsection( 5), error=error) 00214 CALL section_vals_get(rotation_section, explicit=my_subsection( 6), error=error) 00215 CALL section_vals_get(dfunct_section, explicit=my_subsection( 7), error=error) 00216 CALL section_vals_get(qparm_section, explicit=my_subsection( 8), error=error) 00217 CALL section_vals_get(hydronium_section, explicit=my_subsection( 9), error=error) 00218 ! These are just special cases since they are not present in their own defition of COLVARS 00219 IF (ASSOCIATED(reaction_path_section)) THEN 00220 CALL section_vals_get(reaction_path_section,& 00221 explicit=my_subsection(10), error=error) 00222 END IF 00223 IF (ASSOCIATED(distance_from_path_section)) THEN 00224 CALL section_vals_get(distance_from_path_section,& 00225 explicit=my_subsection(16), error=error) 00226 END IF 00227 IF (ASSOCIATED(combine_section)) THEN 00228 CALL section_vals_get(combine_section, explicit=my_subsection(11), error=error) 00229 END IF 00230 CALL section_vals_get(population_section, explicit=my_subsection(12), error=error) 00231 CALL section_vals_get(plane_plane_angle_section,& 00232 explicit=my_subsection(13), error=error) 00233 CALL section_vals_get(gyration_section, explicit=my_subsection(14), error=error) 00234 CALL section_vals_get(rmsd_section, explicit=my_subsection(15), error=error) 00235 CALL section_vals_get(xyz_diag_section, explicit=my_subsection(17), error=error) 00236 CALL section_vals_get(xyz_outerdiag_section,explicit=my_subsection(18), error=error) 00237 CALL section_vals_get(u_section, explicit=my_subsection(19), error=error) 00238 CALL section_vals_get(Wc_section, explicit=my_subsection(20), error=error) 00239 CALL section_vals_get(HBP_section, explicit=my_subsection(21), error=error) 00240 CALL section_vals_get(ring_puckering_section,& 00241 explicit=my_subsection(22), error=error) 00242 CALL section_vals_get(mindist_section, explicit=my_subsection(23), error=error) 00243 00244 ! Only one colvar can be present 00245 CPPostcondition(COUNT(my_subsection) == 1,cp_failure_level,routinep,error,failure) 00246 CPPostcondition(.NOT.ASSOCIATED(colvar),cp_failure_level,routinep,error,failure) 00247 00248 IF (my_subsection(1)) THEN 00249 ! Distance 00250 wrk_section => distance_section 00251 CALL colvar_create(colvar, dist_colvar_id, error) 00252 CALL colvar_check_points(colvar, distance_section, error) 00253 CALL section_vals_val_get(distance_section,"ATOMS",i_vals=iatms,error=error) 00254 colvar%dist_param%i_at = iatms(1) 00255 colvar%dist_param%j_at = iatms(2) 00256 CALL section_vals_val_get(distance_section,"AXIS",i_val=colvar%dist_param%axis_id,error=error) 00257 ELSE IF (my_subsection(2)) THEN 00258 ! Angle 00259 wrk_section => angle_section 00260 CALL colvar_create(colvar, angle_colvar_id, error) 00261 CALL colvar_check_points(colvar, angle_section, error) 00262 CALL section_vals_val_get(angle_section,"ATOMS",i_vals=iatms,error=error) 00263 colvar%angle_param%i_at_angle = iatms 00264 ELSE IF (my_subsection(3)) THEN 00265 ! Torsion 00266 wrk_section => torsion_section 00267 CALL colvar_create(colvar, torsion_colvar_id, error) 00268 CALL colvar_check_points(colvar, torsion_section, error) 00269 CALL section_vals_val_get(torsion_section,"ATOMS",i_vals=iatms,error=error) 00270 colvar%torsion_param%i_at_tors = iatms 00271 colvar%torsion_param%o0 = 0.0_dp 00272 ELSE IF (my_subsection(4)) THEN 00273 ! Coordination 00274 wrk_section => coordination_section 00275 CALL colvar_create(colvar, coord_colvar_id, error) 00276 CALL colvar_check_points(colvar, coordination_section, error) 00277 NULLIFY(colvar%coord_param%i_at_from, colvar%coord_param%c_kinds_from) 00278 NULLIFY(colvar%coord_param%i_at_to, colvar%coord_param%c_kinds_to) 00279 NULLIFY(colvar%coord_param%i_at_to_b, colvar%coord_param%c_kinds_to_b) 00280 ! This section can be repeated 00281 CALL section_vals_val_get(coordination_section,"ATOMS_FROM",n_rep_val=n_var,error=error) 00282 ndim = 0 00283 IF (n_var /= 0) THEN 00284 ! INDEX LIST 00285 DO k = 1, n_var 00286 CALL section_vals_val_get(coordination_section,"ATOMS_FROM",i_rep_val=k,i_vals=iatms,error=error) 00287 CALL reallocate(colvar%coord_param%i_at_from,1, ndim+SIZE(iatms)) 00288 colvar%coord_param%i_at_from(ndim+1:ndim+SIZE(iatms)) = iatms 00289 ndim = ndim + SIZE(iatms) 00290 END DO 00291 colvar%coord_param%n_atoms_from = ndim 00292 colvar%coord_param%use_kinds_from = .FALSE. 00293 ELSE 00294 ! KINDS 00295 CALL section_vals_val_get(coordination_section,"KINDS_FROM",n_rep_val=n_var,error=error) 00296 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 00297 DO k = 1, n_var 00298 CALL section_vals_val_get(coordination_section,"KINDS_FROM",i_rep_val=k,c_vals=c_kinds,error=error) 00299 CALL reallocate(colvar%coord_param%c_kinds_from,1, ndim+SIZE(c_kinds)) 00300 colvar%coord_param%c_kinds_from(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00301 ndim = ndim + SIZE(c_kinds) 00302 END DO 00303 colvar%coord_param%n_atoms_from = 0 00304 colvar%coord_param%use_kinds_from = .TRUE. 00305 ! Uppercase the label 00306 DO k = 1, ndim 00307 CALL uppercase(colvar%coord_param%c_kinds_from(k)) 00308 END DO 00309 END IF 00310 ! This section can be repeated 00311 CALL section_vals_val_get(coordination_section,"ATOMS_TO",n_rep_val=n_var,error=error) 00312 ndim = 0 00313 IF (n_var /= 0) THEN 00314 ! INDEX LIST 00315 DO k = 1, n_var 00316 CALL section_vals_val_get(coordination_section,"ATOMS_TO",i_rep_val=k,i_vals=iatms,error=error) 00317 CALL reallocate(colvar%coord_param%i_at_to,1, ndim+SIZE(iatms)) 00318 colvar%coord_param%i_at_to(ndim+1:ndim+SIZE(iatms)) = iatms 00319 ndim = ndim + SIZE(iatms) 00320 END DO 00321 colvar%coord_param%n_atoms_to = ndim 00322 colvar%coord_param%use_kinds_to = .FALSE. 00323 ELSE 00324 ! KINDS 00325 CALL section_vals_val_get(coordination_section,"KINDS_TO",n_rep_val=n_var,error=error) 00326 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 00327 DO k = 1, n_var 00328 CALL section_vals_val_get(coordination_section,"KINDS_TO",i_rep_val=k,c_vals=c_kinds,error=error) 00329 CALL reallocate(colvar%coord_param%c_kinds_to,1, ndim+SIZE(c_kinds)) 00330 colvar%coord_param%c_kinds_to(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00331 ndim = ndim + SIZE(c_kinds) 00332 END DO 00333 colvar%coord_param%n_atoms_to = 0 00334 colvar%coord_param%use_kinds_to = .TRUE. 00335 ! Uppercase the label 00336 DO k = 1, ndim 00337 CALL uppercase(colvar%coord_param%c_kinds_to(k)) 00338 END DO 00339 END IF 00340 ! Let's finish reading the other parameters 00341 CALL section_vals_val_get(coordination_section,"R0",r_val=colvar%coord_param%r_0,error=error) 00342 CALL section_vals_val_get(coordination_section,"NN",i_val=colvar%coord_param%nncrd,error=error) 00343 CALL section_vals_val_get(coordination_section,"ND",i_val=colvar%coord_param%ndcrd,error=error) 00344 ! This section can be repeated 00345 CALL section_vals_val_get(coordination_section,"ATOMS_TO_B",n_rep_val=n_var,error=error) 00346 CALL section_vals_val_get(coordination_section,"KINDS_TO_B",n_rep_val=n_var_k,error=error) 00347 ndim = 0 00348 IF (n_var /= 0 .OR. n_var_k /= 0) THEN 00349 colvar%coord_param%do_chain = .TRUE. 00350 IF (n_var /= 0) THEN 00351 ! INDEX LIST 00352 DO k = 1, n_var 00353 CALL section_vals_val_get(coordination_section,"ATOMS_TO_B",i_rep_val=k,i_vals=iatms,error=error) 00354 CALL reallocate(colvar%coord_param%i_at_to_b,1, ndim+SIZE(iatms)) 00355 colvar%coord_param%i_at_to_b(ndim+1:ndim+SIZE(iatms)) = iatms 00356 ndim = ndim + SIZE(iatms) 00357 END DO 00358 colvar%coord_param%n_atoms_to_b = ndim 00359 colvar%coord_param%use_kinds_to_b = .FALSE. 00360 ELSE 00361 ! KINDS 00362 CALL section_vals_val_get(coordination_section,"KINDS_TO_B",n_rep_val=n_var_k,error=error) 00363 CPPostcondition(n_var_k>0,cp_failure_level,routinep,error,failure) 00364 DO k = 1, n_var_k 00365 CALL section_vals_val_get(coordination_section,"KINDS_TO_B",i_rep_val=k,c_vals=c_kinds,error=error) 00366 CALL reallocate(colvar%coord_param%c_kinds_to_b,1, ndim+SIZE(c_kinds)) 00367 colvar%coord_param%c_kinds_to_b(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00368 ndim = ndim + SIZE(c_kinds) 00369 END DO 00370 colvar%coord_param%n_atoms_to_b = 0 00371 colvar%coord_param%use_kinds_to_b = .TRUE. 00372 ! Uppercase the label 00373 DO k = 1, ndim 00374 CALL uppercase(colvar%coord_param%c_kinds_to_b(k)) 00375 END DO 00376 END IF 00377 ! Let's finish reading the other parameters 00378 CALL section_vals_val_get(coordination_section,"R0_B",r_val=colvar%coord_param%r_0_b,error=error) 00379 CALL section_vals_val_get(coordination_section,"NN_B",i_val=colvar%coord_param%nncrd_b,error=error) 00380 CALL section_vals_val_get(coordination_section,"ND_B",i_val=colvar%coord_param%ndcrd_b,error=error) 00381 ELSE 00382 colvar%coord_param%do_chain = .FALSE. 00383 colvar%coord_param%n_atoms_to_b = 0 00384 colvar%coord_param%use_kinds_to_b = .FALSE. 00385 NULLIFY(colvar%coord_param%i_at_to_b) 00386 NULLIFY(colvar%coord_param%c_kinds_to_b) 00387 colvar%coord_param%nncrd_b = 0 00388 colvar%coord_param%ndcrd_b = 0 00389 colvar%coord_param%r_0_b = 0._dp 00390 END IF 00391 00392 ELSE IF (my_subsection(5)) THEN 00393 ! Distance point from plane 00394 wrk_section => plane_dist_section 00395 CALL colvar_create(colvar, plane_distance_colvar_id, error) 00396 CALL colvar_check_points(colvar, plane_dist_section, error) 00397 CALL section_vals_val_get(plane_dist_section,"ATOMS_PLANE",i_vals=iatms,error=error) 00398 CPPostcondition(SIZE(iatms) == 3,cp_failure_level,routinep,error,failure) 00399 colvar%plane_distance_param%plane = iatms 00400 CALL section_vals_val_get(plane_dist_section,"ATOM_POINT",i_val=iatm,error=error) 00401 colvar%plane_distance_param%point = iatm 00402 CALL section_vals_val_get(plane_dist_section,"PBC",l_val=colvar%plane_distance_param%use_pbc,error=error) 00403 ELSE IF (my_subsection(6)) THEN 00404 ! Rotation colvar of a segment w.r.t. another segment 00405 wrk_section => rotation_section 00406 CALL colvar_create(colvar, rotation_colvar_id, error) 00407 CALL colvar_check_points(colvar, rotation_section, error) 00408 CALL section_vals_val_get(rotation_section,"P1_BOND1",i_val=colvar%rotation_param%i_at1_bond1,error=error) 00409 CALL section_vals_val_get(rotation_section,"P2_BOND1",i_val=colvar%rotation_param%i_at2_bond1,error=error) 00410 CALL section_vals_val_get(rotation_section,"P1_BOND2",i_val=colvar%rotation_param%i_at1_bond2,error=error) 00411 CALL section_vals_val_get(rotation_section,"P2_BOND2",i_val=colvar%rotation_param%i_at2_bond2,error=error) 00412 ELSE IF (my_subsection(7)) THEN 00413 ! Difference of two distances 00414 wrk_section => dfunct_section 00415 CALL colvar_create(colvar, dfunct_colvar_id, error) 00416 CALL colvar_check_points(colvar, dfunct_section, error) 00417 CALL section_vals_val_get(dfunct_section,"ATOMS",i_vals=iatms,error=error) 00418 colvar%dfunct_param%i_at_dfunct = iatms 00419 CALL section_vals_val_get(dfunct_section,"COEFFICIENT",r_val=colvar%dfunct_param%coeff,error=error) 00420 CALL section_vals_val_get(dfunct_section,"PBC",l_val=colvar%dfunct_param%use_pbc,error=error) 00421 ELSE IF (my_subsection(8)) THEN 00422 ! Q Parameter 00423 wrk_section => qparm_section 00424 CALL colvar_create(colvar, qparm_colvar_id, error) 00425 CALL colvar_check_points(colvar, qparm_section, error) 00426 CALL section_vals_val_get(qparm_section,"RCUT",r_val=colvar%qparm_param%rcut,error=error) 00427 CALL section_vals_val_get(qparm_section,"ALPHA",r_val=colvar%qparm_param%alpha,error=error) 00428 CALL section_vals_val_get(qparm_section,"L",i_val=colvar%qparm_param%l,error=error) 00429 NULLIFY(colvar%qparm_param%i_at_from) 00430 NULLIFY(colvar%qparm_param%i_at_to) 00431 CALL section_vals_val_get(qparm_section,"ATOMS_FROM",n_rep_val=n_var,error=error) 00432 ndim = 0 00433 DO k = 1, n_var 00434 CALL section_vals_val_get(qparm_section,"ATOMS_FROM",i_rep_val=k,i_vals=iatms,error=error) 00435 CALL reallocate(colvar%qparm_param%i_at_from,1, ndim+SIZE(iatms)) 00436 colvar%qparm_param%i_at_from(ndim+1:ndim+SIZE(iatms)) = iatms 00437 ndim = ndim + SIZE(iatms) 00438 END DO 00439 colvar%qparm_param%n_atoms_from = ndim 00440 ! This section can be repeated 00441 CALL section_vals_val_get(qparm_section,"ATOMS_TO",n_rep_val=n_var,error=error) 00442 ndim = 0 00443 DO k = 1, n_var 00444 CALL section_vals_val_get(qparm_section,"ATOMS_TO",i_rep_val=k,i_vals=iatms,error=error) 00445 CALL reallocate(colvar%qparm_param%i_at_to,1, ndim+SIZE(iatms)) 00446 colvar%qparm_param%i_at_to(ndim+1:ndim+SIZE(iatms)) = iatms 00447 ndim = ndim + SIZE(iatms) 00448 END DO 00449 colvar%qparm_param%n_atoms_to = ndim 00450 ELSE IF (my_subsection(9)) THEN 00451 ! Hydronium 00452 CALL colvar_create(colvar,hydronium_colvar_id, error) 00453 NULLIFY(colvar%hydronium_param%i_oxygens) 00454 NULLIFY(colvar%hydronium_param%i_hydrogens) 00455 00456 CALL section_vals_val_get(hydronium_section,"OXYGENS",n_rep_val=n_var,error=error) 00457 ndim = 0 00458 DO k = 1, n_var 00459 CALL section_vals_val_get(hydronium_section,"OXYGENS",i_vals=iatms,error=error) 00460 CALL reallocate(colvar%hydronium_param%i_oxygens,1,ndim+SIZE(iatms)) 00461 colvar%hydronium_param%i_oxygens(ndim+1:ndim+SIZE(iatms)) = iatms 00462 ndim = ndim + SIZE(iatms) 00463 END DO 00464 colvar%hydronium_param%n_oxygens = ndim 00465 00466 CALL section_vals_val_get(hydronium_section,"HYDROGENS",n_rep_val=n_var,error=error) 00467 ndim = 0 00468 DO k = 1, n_var 00469 CALL section_vals_val_get(hydronium_section,"HYDROGENS",i_vals=iatms,error=error) 00470 CALL reallocate(colvar%hydronium_param%i_hydrogens,1,ndim+SIZE(iatms)) 00471 colvar%hydronium_param%i_hydrogens(ndim+1:ndim+SIZE(iatms)) = iatms 00472 ndim = ndim + SIZE(iatms) 00473 END DO 00474 colvar%hydronium_param%n_hydrogens = ndim 00475 00476 CALL section_vals_val_get(hydronium_section,"ROO",r_val=colvar%hydronium_param%r_OO,error=error) 00477 CALL section_vals_val_get(hydronium_section,"ROH",r_val=colvar%hydronium_param%r_OH,error=error) 00478 CALL section_vals_val_get(hydronium_section,"pNH",i_val=colvar%hydronium_param%pnh,error=error) 00479 CALL section_vals_val_get(hydronium_section,"qNH",i_val=colvar%hydronium_param%qnh,error=error) 00480 CALL section_vals_val_get(hydronium_section,"pNO",i_val=colvar%hydronium_param%pno,error=error) 00481 CALL section_vals_val_get(hydronium_section,"qNO",i_val=colvar%hydronium_param%qno,error=error) 00482 CALL section_vals_val_get(hydronium_section,"p",i_val=colvar%hydronium_param%p,error=error) 00483 CALL section_vals_val_get(hydronium_section,"q",i_val=colvar%hydronium_param%q,error=error) 00484 CALL section_vals_val_get(hydronium_section,"NH",r_val=colvar%hydronium_param%nh,error=error) 00485 CALL section_vals_val_get(hydronium_section,"LAMBDA",r_val=colvar%hydronium_param%lambda,error=error) 00486 ELSE IF(my_subsection(10) .OR. my_subsection(16))THEN 00487 !reaction path or distance from reaction path 00488 IF (my_subsection(10) ) THEN 00489 path_section => reaction_path_section 00490 CALL colvar_create(colvar,reaction_path_colvar_id, error) 00491 fmid="POS" 00492 ifunc=1 00493 ELSE IF (my_subsection(16) ) THEN 00494 path_section => distance_from_path_section 00495 CALL colvar_create(colvar,distance_from_path_colvar_id, error) 00496 fmid="DIS" 00497 ifunc=2 00498 END IF 00499 colvar%use_points=.FALSE. 00500 CALL section_vals_val_get(path_section,"LAMBDA",r_val=colvar%reaction_path_param%lambda,error=error) 00501 CALL section_vals_val_get(path_section,"DISTANCES_RMSD",l_val=colvar%reaction_path_param%dist_rmsd,& 00502 error=error) 00503 CALL section_vals_val_get(path_section,"RMSD",l_val=colvar%reaction_path_param%rmsd,& 00504 error=error) 00505 IF(colvar%reaction_path_param%dist_rmsd .AND. colvar%reaction_path_param%rmsd)THEN 00506 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00507 "CV REACTION PATH: only one between DISTANCES_RMSD and RMSD can be used "//& 00508 CPSourceFileRef) 00509 END IF 00510 IF(colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd) THEN 00511 NULLIFY(colvar%reaction_path_param%i_rmsd, colvar%reaction_path_param%r_ref) 00512 frame_section => section_vals_get_subs_vals(path_section,"FRAME",error=error) 00513 CALL section_vals_get(frame_section,n_repetition=nr_frame,error=error) 00514 00515 colvar%reaction_path_param%nr_frames=nr_frame 00516 CALL read_frames(frame_section,para_env,nr_frame,colvar%reaction_path_param%r_ref,& 00517 colvar%reaction_path_param%n_components,error=error) 00518 CALL section_vals_val_get(path_section,"SUBSET_TYPE",i_val=colvar%reaction_path_param%subset,& 00519 error=error) 00520 IF (colvar%reaction_path_param%subset==rmsd_all) THEN 00521 ALLOCATE(colvar%reaction_path_param%i_rmsd(colvar%reaction_path_param%n_components),STAT=stat) 00522 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00523 DO i = 1, colvar%reaction_path_param%n_components 00524 colvar%reaction_path_param%i_rmsd(i) = i 00525 END DO 00526 ELSE IF (colvar%reaction_path_param%subset==rmsd_list) THEN 00527 ! This section can be repeated 00528 CALL section_vals_val_get(path_section,"ATOMS",n_rep_val=n_var,error=error) 00529 ndim = 0 00530 IF (n_var /= 0) THEN 00531 ! INDEX LIST 00532 DO k = 1, n_var 00533 CALL section_vals_val_get(path_section,"ATOMS",i_rep_val=k,i_vals=iatms,error=error) 00534 CALL reallocate(colvar%reaction_path_param%i_rmsd,1, ndim+SIZE(iatms)) 00535 colvar%reaction_path_param%i_rmsd(ndim+1:ndim+SIZE(iatms)) = iatms 00536 ndim = ndim + SIZE(iatms) 00537 END DO 00538 colvar%reaction_path_param%n_components = ndim 00539 ELSE 00540 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00541 "CV REACTION PATH: if SUBSET_TYPE=LIST a list of atoms needs to be provided "//& 00542 CPSourceFileRef) 00543 END IF 00544 END IF 00545 00546 CALL section_vals_val_get(path_section,"ALIGN_FRAMES",l_val=colvar%reaction_path_param%align_frames,& 00547 error=error) 00548 ELSE 00549 colvar_subsection => section_vals_get_subs_vals(path_section,"COLVAR",error=error) 00550 CALL section_vals_get(colvar_subsection,n_repetition=ncol,error=error) 00551 ALLOCATE(colvar%reaction_path_param%colvar_p(ncol),stat=stat) 00552 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00553 IF (ncol>0) THEN 00554 DO i= 1, ncol 00555 NULLIFY(colvar%reaction_path_param%colvar_p(i)%colvar) 00556 CALL colvar_read(colvar%reaction_path_param%colvar_p(i)%colvar,i,colvar_subsection, para_env,& 00557 error=error) 00558 ENDDO 00559 ELSE 00560 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00561 "CV REACTION PATH: the number of CV to define the path must be >0 "//& 00562 CPSourceFileRef) 00563 ENDIF 00564 colvar%reaction_path_param%n_components=ncol 00565 NULLIFY(range) 00566 CALL section_vals_val_get(path_section,"RANGE",r_vals=range,error=error) 00567 CALL section_vals_val_get(path_section,"STEP_SIZE",r_val=colvar%reaction_path_param%step_size,error=error) 00568 iend=CEILING(MAX(RANGE(1),RANGE(2))/colvar%reaction_path_param%step_size) 00569 istart=FLOOR(MIN(RANGE(1),RANGE(2))/colvar%reaction_path_param%step_size) 00570 colvar%reaction_path_param%function_bounds(1)=istart 00571 colvar%reaction_path_param%function_bounds(2)=iend 00572 colvar%reaction_path_param%nr_frames= 2 !iend - istart + 1 00573 ALLOCATE(colvar%reaction_path_param%f_vals(ncol,istart:iend),stat=stat) 00574 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00575 CALL section_vals_val_get(path_section,"VARIABLE",c_vals=my_par,i_rep_val=1,error=error) 00576 CALL section_vals_val_get(path_section,"FUNCTION",n_rep_val=ncol,error=error) 00577 check = (ncol==SIZE(colvar%reaction_path_param%colvar_p)) 00578 CPPostcondition(check,cp_failure_level,routinep,error,failure) 00579 CALL initf(ncol) 00580 DO i=1,ncol 00581 CALL section_vals_val_get(path_section,"FUNCTION",c_val=path_function,i_rep_val=i,error=error) 00582 CALL compress(path_function, full=.TRUE.) 00583 CALL parsef(i,TRIM(path_function),my_par) 00584 DO j=istart,iend 00585 my_val=REAL(j,kind=dp)*colvar%reaction_path_param%step_size 00586 colvar%reaction_path_param%f_vals(i,j)=evalf(i,my_val) 00587 END DO 00588 END DO 00589 CALL finalizef() 00590 00591 iw1= cp_print_key_unit_nr(logger,path_section,& 00592 "MAP",middle_name=fmid,extension=".dat",file_status="REPLACE",error=error) 00593 IF(iw1>0)THEN 00594 CALL section_vals_val_get(path_section,"MAP%GRID_SPACING",n_rep_val=ncol,error=error) 00595 ALLOCATE(grid_sp(ncol),stat=stat) 00596 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00597 DO i=1,ncol 00598 CALL section_vals_val_get(path_section,"MAP%GRID_SPACING",r_val=grid_sp(i),error=error) 00599 END DO 00600 CALL section_vals_val_get(path_section,"MAP%RANGE",n_rep_val=ncol,error=error) 00601 CPPostcondition(ncol ==SIZE(grid_sp) ,cp_failure_level,routinep,error,failure) 00602 ALLOCATE(p_range(2,ncol),stat=stat) 00603 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00604 ALLOCATE(p_bounds(2,ncol),stat=stat) 00605 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00606 DO i=1,ncol 00607 CALL section_vals_val_get(path_section,"MAP%RANGE",r_vals=g_range,error=error) 00608 p_range(:,i)=g_range(:) 00609 p_bounds(2,i)=CEILING( MAX(p_range(1,i),p_range(2,i))/grid_sp(i)) 00610 p_bounds(1,i)=FLOOR(MIN(p_range(1,i) , p_range(2,i))/grid_sp(i)) 00611 END DO 00612 ALLOCATE(s1v(2,istart:iend),stat=stat) 00613 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00614 ALLOCATE(s1(2),stat=stat) 00615 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00616 ALLOCATE(grid_point(ncol),stat=stat) 00617 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00618 v_count=0 00619 kk=rec_eval_grid(iw1,ncol,colvar%reaction_path_param%f_vals,v_count,& 00620 grid_point,grid_sp,colvar%reaction_path_param%step_size,istart,& 00621 iend,s1v,s1,p_bounds,colvar%reaction_path_param%lambda,ifunc=ifunc,& 00622 nconf=colvar%reaction_path_param%nr_frames) 00623 DEALLOCATE(grid_sp,stat=stat) 00624 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00625 DEALLOCATE(p_range,stat=stat) 00626 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00627 DEALLOCATE(p_bounds,stat=stat) 00628 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00629 DEALLOCATE(s1v,stat=stat) 00630 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00631 DEALLOCATE(s1,stat=stat) 00632 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00633 DEALLOCATE(grid_point,stat=stat) 00634 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00635 END IF 00636 CALL cp_print_key_finished_output(iw1,logger,path_section,& 00637 "MAP", error=error) 00638 END IF 00639 00640 ELSE IF(my_subsection(11))THEN 00641 ! combine colvar 00642 CALL colvar_create(colvar,combine_colvar_id, error) 00643 colvar%use_points=.FALSE. 00644 colvar_subsection => section_vals_get_subs_vals(combine_section,"COLVAR",error=error) 00645 CALL section_vals_get(colvar_subsection,n_repetition=ncol,error=error) 00646 ALLOCATE(colvar%combine_cvs_param%colvar_p(ncol),stat=stat) 00647 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00648 ! In case we need to print some information.. 00649 iw = cp_print_key_unit_nr(logger,colvar_section,& 00650 "PRINT%PROGRAM_RUN_INFO",extension=".colvarLog",error=error) 00651 IF (iw>0) THEN 00652 WRITE ( iw, '( A )')' '//& 00653 '**********************************************************************' 00654 WRITE ( iw, '( A,I8)' ) ' COLVARS| COLVAR INPUT INDEX: ',icol 00655 WRITE ( iw, '( A,T49,4I8)' ) ' COLVARS| COMBINATION OF THE FOLOWING COLVARS:' 00656 END IF 00657 CALL cp_print_key_finished_output(iw,logger,colvar_section,& 00658 "PRINT%PROGRAM_RUN_INFO", error=error) 00659 ! Parsing the real COLVARs 00660 DO i= 1, ncol 00661 NULLIFY(colvar%combine_cvs_param%colvar_p(i)%colvar) 00662 CALL colvar_read(colvar%combine_cvs_param%colvar_p(i)%colvar,i,colvar_subsection, para_env, error=error) 00663 END DO 00664 ! Function definition 00665 CALL section_vals_val_get(combine_section,"FUNCTION",c_val=colvar%combine_cvs_param%function,error=error) 00666 CALL compress(colvar%combine_cvs_param%function, full=.TRUE.) 00667 ! Variables 00668 CALL section_vals_val_get(combine_section,"VARIABLES",c_vals=my_par,error=error) 00669 ALLOCATE(colvar%combine_cvs_param%variables(SIZE(my_par)),stat=stat) 00670 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00671 colvar%combine_cvs_param%variables = my_par 00672 ! Check that the number of COLVAR provided is equal to the number of variables.. 00673 CALL cp_assert(SIZE(my_par)==ncol,cp_fatal_level,cp_assertion_failed,routineP,& 00674 "Number of defined COLVAR for COMBINE_COLVAR is different from the "//& 00675 "number of variables! It is not possible to define COLVARs in a COMBINE_COLVAR "//& 00676 "and avoid their usage in the combininig function!"//& 00677 CPSourceFileRef) 00678 ! Parameters 00679 ALLOCATE(colvar%combine_cvs_param%c_parameters(0),stat=stat) 00680 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00681 CALL section_vals_val_get(combine_section,"PARAMETERS",n_rep_val=ncol,error=error) 00682 DO i = 1,ncol 00683 isize = SIZE(colvar%combine_cvs_param%c_parameters) 00684 CALL section_vals_val_get(combine_section,"PARAMETERS",c_vals=my_par,i_rep_val=i,error=error) 00685 CALL reallocate(colvar%combine_cvs_param%c_parameters,1,isize+SIZE(my_par)) 00686 colvar%combine_cvs_param%c_parameters(isize+1:isize+SIZE(my_par)) = my_par 00687 END DO 00688 ALLOCATE(colvar%combine_cvs_param%v_parameters(0),stat=stat) 00689 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00690 CALL section_vals_val_get(combine_section,"VALUES",n_rep_val=ncol,error=error) 00691 DO i = 1,ncol 00692 isize = SIZE(colvar%combine_cvs_param%v_parameters) 00693 CALL section_vals_val_get(combine_section,"VALUES",r_vals=my_vals,i_rep_val=i,error=error) 00694 CALL reallocate(colvar%combine_cvs_param%v_parameters,1,isize+SIZE(my_vals)) 00695 colvar%combine_cvs_param%v_parameters(isize+1:isize+SIZE(my_vals)) = my_vals 00696 END DO 00697 ! Info on derivative evaluation 00698 CALL section_vals_val_get(combine_section,"DX",r_val=colvar%combine_cvs_param%dx,error=error) 00699 CALL section_vals_val_get(combine_section,"ERROR_LIMIT",r_val=colvar%combine_cvs_param%lerr,& 00700 error=error) 00701 ELSE IF (my_subsection(12)) THEN 00702 ! Population 00703 wrk_section => population_section 00704 CALL colvar_create(colvar, population_colvar_id, error) 00705 CALL colvar_check_points(colvar, population_section, error) 00706 00707 NULLIFY(colvar%population_param%i_at_from, colvar%population_param%c_kinds_from) 00708 NULLIFY(colvar%population_param%i_at_to, colvar%population_param%c_kinds_to) 00709 ! This section can be repeated 00710 00711 CALL section_vals_val_get(population_section,"ATOMS_FROM",n_rep_val=n_var,error=error) 00712 ndim = 0 00713 IF (n_var /= 0) THEN 00714 ! INDEX LIST 00715 DO k = 1, n_var 00716 CALL section_vals_val_get(population_section,"ATOMS_FROM",i_rep_val=k,i_vals=iatms,error=error) 00717 CALL reallocate(colvar%population_param%i_at_from,1, ndim+SIZE(iatms)) 00718 colvar%population_param%i_at_from(ndim+1:ndim+SIZE(iatms)) = iatms 00719 ndim = ndim + SIZE(iatms) 00720 END DO 00721 colvar%population_param%n_atoms_from = ndim 00722 colvar%population_param%use_kinds_from = .FALSE. 00723 ELSE 00724 ! KINDS 00725 CALL section_vals_val_get(population_section,"KINDS_FROM",n_rep_val=n_var,error=error) 00726 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 00727 DO k = 1, n_var 00728 CALL section_vals_val_get(population_section,"KINDS_FROM",i_rep_val=k,c_vals=c_kinds,error=error) 00729 CALL reallocate(colvar%population_param%c_kinds_from,1, ndim+SIZE(c_kinds)) 00730 colvar%population_param%c_kinds_from(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00731 ndim = ndim + SIZE(c_kinds) 00732 END DO 00733 colvar%population_param%n_atoms_from = 0 00734 colvar%population_param%use_kinds_from = .TRUE. 00735 ! Uppercase the label 00736 DO k = 1, ndim 00737 CALL uppercase(colvar%population_param%c_kinds_from(k)) 00738 END DO 00739 END IF 00740 ! This section can be repeated 00741 CALL section_vals_val_get(population_section,"ATOMS_TO",n_rep_val=n_var,error=error) 00742 ndim = 0 00743 IF (n_var /= 0) THEN 00744 ! INDEX LIST 00745 DO k = 1, n_var 00746 CALL section_vals_val_get(population_section,"ATOMS_TO",i_rep_val=k,i_vals=iatms,error=error) 00747 CALL reallocate(colvar%population_param%i_at_to,1, ndim+SIZE(iatms)) 00748 colvar%population_param%i_at_to(ndim+1:ndim+SIZE(iatms)) = iatms 00749 ndim = ndim + SIZE(iatms) 00750 END DO 00751 colvar%population_param%n_atoms_to = ndim 00752 colvar%population_param%use_kinds_to = .FALSE. 00753 ELSE 00754 ! KINDS 00755 CALL section_vals_val_get(population_section,"KINDS_TO",n_rep_val=n_var,error=error) 00756 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 00757 DO k = 1, n_var 00758 CALL section_vals_val_get(population_section,"KINDS_TO",i_rep_val=k,c_vals=c_kinds,error=error) 00759 CALL reallocate(colvar%population_param%c_kinds_to,1, ndim+SIZE(c_kinds)) 00760 colvar%population_param%c_kinds_to(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00761 ndim = ndim + SIZE(c_kinds) 00762 END DO 00763 colvar%population_param%n_atoms_to = 0 00764 colvar%population_param%use_kinds_to = .TRUE. 00765 ! Uppercase the label 00766 DO k = 1, ndim 00767 CALL uppercase(colvar%population_param%c_kinds_to(k)) 00768 END DO 00769 END IF 00770 ! Let's finish reading the other parameters 00771 CALL section_vals_val_get(population_section,"R0",r_val=colvar%population_param%r_0,error=error) 00772 CALL section_vals_val_get(population_section,"NN",i_val=colvar%population_param%nncrd,error=error) 00773 CALL section_vals_val_get(population_section,"ND",i_val=colvar%population_param%ndcrd,error=error) 00774 CALL section_vals_val_get(population_section,"N0",i_val=colvar%population_param%n0,error=error) 00775 CALL section_vals_val_get(population_section,"SIGMA",r_val=colvar%population_param%sigma,error=error) 00776 ELSE IF (my_subsection(13)) THEN 00777 ! Angle between two planes 00778 wrk_section => plane_plane_angle_section 00779 CALL colvar_create(colvar, plane_plane_angle_colvar_id, error) 00780 CALL colvar_check_points(colvar, plane_plane_angle_section, error) 00781 ! Read the specification of the two planes 00782 plane_sections => section_vals_get_subs_vals(plane_plane_angle_section,"PLANE",error=error) 00783 CALL section_vals_get(plane_sections, n_repetition=n_var, error=error) 00784 CALL cp_assert(n_var==2,cp_fatal_level,cp_assertion_failed,routineP,& 00785 "PLANE_PLANE_ANGLE Colvar section: Two PLANE sections must be provided!"//& 00786 CPSourceFileRef) 00787 ! Plane 1 00788 CALL section_vals_val_get(plane_sections,"DEF_TYPE",i_rep_section=1,& 00789 i_val=colvar%plane_plane_angle_param%plane1%type_of_def,error=error) 00790 IF (colvar%plane_plane_angle_param%plane1%type_of_def==plane_def_vec) THEN 00791 CALL section_vals_val_get(plane_sections,"NORMAL_VECTOR",i_rep_section=1,& 00792 r_vals=s1,error=error) 00793 colvar%plane_plane_angle_param%plane1%normal_vec = s1 00794 ELSE 00795 CALL section_vals_val_get(plane_sections,"ATOMS",i_rep_section=1,& 00796 i_vals=iatms,error=error) 00797 colvar%plane_plane_angle_param%plane1%points = iatms 00798 END IF 00799 00800 ! Plane 2 00801 CALL section_vals_val_get(plane_sections,"DEF_TYPE",i_rep_section=2,& 00802 i_val=colvar%plane_plane_angle_param%plane2%type_of_def,error=error) 00803 IF (colvar%plane_plane_angle_param%plane2%type_of_def==plane_def_vec) THEN 00804 CALL section_vals_val_get(plane_sections,"NORMAL_VECTOR",i_rep_section=2,& 00805 r_vals=s1,error=error) 00806 colvar%plane_plane_angle_param%plane2%normal_vec = s1 00807 ELSE 00808 CALL section_vals_val_get(plane_sections,"ATOMS",i_rep_section=2,& 00809 i_vals=iatms,error=error) 00810 colvar%plane_plane_angle_param%plane2%points = iatms 00811 END IF 00812 ELSE IF (my_subsection(14)) THEN 00813 ! Gyration Radius 00814 wrk_section => gyration_section 00815 CALL colvar_create(colvar, gyration_colvar_id, error) 00816 CALL colvar_check_points(colvar, gyration_section, error) 00817 00818 NULLIFY(colvar%gyration_param%i_at, colvar%gyration_param%c_kinds) 00819 00820 ! This section can be repeated 00821 CALL section_vals_val_get(gyration_section,"ATOMS",n_rep_val=n_var,error=error) 00822 ndim = 0 00823 IF (n_var /= 0) THEN 00824 ! INDEX LIST 00825 DO k = 1, n_var 00826 CALL section_vals_val_get(gyration_section,"ATOMS",i_rep_val=k,i_vals=iatms,error=error) 00827 CALL reallocate(colvar%gyration_param%i_at,1, ndim+SIZE(iatms)) 00828 colvar%gyration_param%i_at(ndim+1:ndim+SIZE(iatms)) = iatms 00829 ndim = ndim + SIZE(iatms) 00830 END DO 00831 colvar%gyration_param%n_atoms = ndim 00832 colvar%gyration_param%use_kinds = .FALSE. 00833 ELSE 00834 ! KINDS 00835 CALL section_vals_val_get(gyration_section,"KINDS",n_rep_val=n_var,error=error) 00836 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 00837 DO k = 1, n_var 00838 CALL section_vals_val_get(gyration_section,"KINDS",i_rep_val=k,c_vals=c_kinds,error=error) 00839 CALL reallocate(colvar%gyration_param%c_kinds,1, ndim+SIZE(c_kinds)) 00840 colvar%gyration_param%c_kinds(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 00841 ndim = ndim + SIZE(c_kinds) 00842 END DO 00843 colvar%gyration_param%n_atoms = 0 00844 colvar%gyration_param%use_kinds = .TRUE. 00845 ! Uppercase the label 00846 DO k = 1, ndim 00847 CALL uppercase(colvar%gyration_param%c_kinds(k)) 00848 END DO 00849 END IF 00850 ELSE IF (my_subsection(15)) THEN 00851 ! RMSD_AB 00852 wrk_section => rmsd_section 00853 CALL colvar_create(colvar, rmsd_colvar_id, error) 00854 00855 NULLIFY(colvar%rmsd_param%i_rmsd, colvar%rmsd_param%r_ref, colvar%rmsd_param%weights) 00856 00857 frame_section => section_vals_get_subs_vals(rmsd_section,"FRAME",error=error) 00858 CALL section_vals_get(frame_section,n_repetition=nr_frame,error=error) 00859 00860 colvar%rmsd_param%nr_frames=nr_frame 00861 ! Calculation is aborted if reference frame are less than 2 00862 CPPostcondition((nr_frame>=1.AND.nr_frame<=2),cp_failure_level,routineP,error,failure) 00863 CALL read_frames(frame_section,para_env,nr_frame,colvar%rmsd_param%r_ref,& 00864 colvar%rmsd_param%n_atoms,error=error) 00865 00866 ALLOCATE(colvar%rmsd_param%weights(colvar%rmsd_param%n_atoms), STAT=stat) 00867 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00868 colvar%rmsd_param%weights = 0.0_dp 00869 00870 CALL section_vals_val_get(rmsd_section,"SUBSET_TYPE",i_val=colvar%rmsd_param%subset,error=error) 00871 IF (colvar%rmsd_param%subset==rmsd_all) THEN 00872 ALLOCATE(colvar%rmsd_param%i_rmsd(colvar%rmsd_param%n_atoms),STAT=stat) 00873 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00874 DO i = 1, colvar%rmsd_param%n_atoms 00875 colvar%rmsd_param%i_rmsd(i) = i 00876 END DO 00877 ELSE IF (colvar%rmsd_param%subset==rmsd_list) THEN 00878 ! This section can be repeated 00879 CALL section_vals_val_get(rmsd_section,"ATOMS",n_rep_val=n_var,error=error) 00880 ndim = 0 00881 IF (n_var /= 0) THEN 00882 ! INDEX LIST 00883 DO k = 1, n_var 00884 CALL section_vals_val_get(rmsd_section,"ATOMS",i_rep_val=k,i_vals=iatms,error=error) 00885 CALL reallocate(colvar%rmsd_param%i_rmsd,1, ndim+SIZE(iatms)) 00886 colvar%rmsd_param%i_rmsd(ndim+1:ndim+SIZE(iatms)) = iatms 00887 ndim = ndim + SIZE(iatms) 00888 END DO 00889 colvar%rmsd_param%n_atoms = ndim 00890 ELSE 00891 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00892 "CV RMSD: if SUBSET_TYPE=LIST a list of atoms needs to be provided "//& 00893 CPSourceFileRef) 00894 END IF 00895 ELSE IF (colvar%rmsd_param%subset==rmsd_weightlist) THEN 00896 CALL section_vals_val_get(rmsd_section,"ATOMS",n_rep_val=n_var,error=error) 00897 ndim = 0 00898 IF (n_var /= 0) THEN 00899 ! INDEX LIST 00900 DO k = 1, n_var 00901 CALL section_vals_val_get(rmsd_section,"ATOMS",i_rep_val=k,i_vals=iatms,error=error) 00902 CALL reallocate(colvar%rmsd_param%i_rmsd,1, ndim+SIZE(iatms)) 00903 colvar%rmsd_param%i_rmsd(ndim+1:ndim+SIZE(iatms)) = iatms 00904 ndim = ndim + SIZE(iatms) 00905 END DO 00906 colvar%rmsd_param%n_atoms = ndim 00907 ELSE 00908 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00909 "CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of atoms needs to be provided "//& 00910 CPSourceFileRef) 00911 END IF 00912 CALL section_vals_val_get(rmsd_section,"WEIGHTS",n_rep_val=n_var,error=error) 00913 ndim = 0 00914 IF (n_var /= 0) THEN 00915 ! INDEX LIST 00916 DO k = 1, n_var 00917 CALL section_vals_val_get(rmsd_section,"WEIGHTS",i_rep_val=k,r_vals=wei,error=error) 00918 CALL reallocate(weights,1, ndim+SIZE(wei)) 00919 weights(ndim+1:ndim+SIZE(wei)) = wei 00920 ndim = ndim + SIZE(wei) 00921 END DO 00922 CALL cp_assert(ndim==colvar%rmsd_param%n_atoms,cp_fatal_level,cp_assertion_failed,routineP,& 00923 "CV RMSD: list of atoms and list of weights need to contain same number of entries. "//& 00924 CPSourceFileRef) 00925 DO i = 1, ndim 00926 ii = colvar%rmsd_param%i_rmsd(i) 00927 colvar%rmsd_param%weights(ii) = weights(i) 00928 END DO 00929 DEALLOCATE (weights, STAT=stat) 00930 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00931 ELSE 00932 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00933 "CV RMSD: if SUBSET_TYPE=WEIGHT_LIST a list of weights need to be provided. "//& 00934 CPSourceFileRef) 00935 END IF 00936 00937 ELSE 00938 CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,routineP,& 00939 "CV RMSD: unknown SUBSET_TYPE."//& 00940 CPSourceFileRef) 00941 END IF 00942 00943 CALL section_vals_val_get(rmsd_section,"ALIGN_FRAMES",l_val=colvar%rmsd_param%align_frames,& 00944 error=error) 00945 ELSE IF (my_subsection(17)) THEN 00946 ! Work on XYZ positions of atoms 00947 wrk_section => xyz_diag_section 00948 CALL colvar_create(colvar, xyz_diag_colvar_id, error) 00949 CALL colvar_check_points(colvar, wrk_section, error) 00950 CALL section_vals_val_get(wrk_section,"ATOM",i_val=iatm,error=error) 00951 CALL section_vals_val_get(wrk_section,"COMPONENT",i_val=icomponent,error=error) 00952 CALL section_vals_val_get(wrk_section,"PBC",l_val=colvar%xyz_diag_param%use_pbc,error=error) 00953 colvar%xyz_diag_param%i_atom = iatm 00954 colvar%xyz_diag_param%component = icomponent 00955 ELSE IF (my_subsection(18)) THEN 00956 ! Work on the outer diagonal (two atoms A,B) XYZ positions 00957 wrk_section => xyz_outerdiag_section 00958 CALL colvar_create(colvar, xyz_outerdiag_colvar_id, error) 00959 CALL colvar_check_points(colvar, wrk_section, error) 00960 CALL section_vals_val_get(wrk_section,"ATOMS",i_vals=iatms,error=error) 00961 colvar%xyz_outerdiag_param%i_atoms = iatms 00962 CALL section_vals_val_get(wrk_section,"COMPONENT_A",i_val=icomponent,error=error) 00963 colvar%xyz_outerdiag_param%components(1) = icomponent 00964 CALL section_vals_val_get(wrk_section,"COMPONENT_B",i_val=icomponent,error=error) 00965 colvar%xyz_outerdiag_param%components(2) = icomponent 00966 CALL section_vals_val_get(wrk_section,"PBC",l_val=colvar%xyz_outerdiag_param%use_pbc,error=error) 00967 ELSE IF (my_subsection(19)) THEN 00968 ! Energy 00969 wrk_section => u_section 00970 CALL colvar_create(colvar, u_colvar_id, error) 00971 colvar%u_param%mixed_energy_section => section_vals_get_subs_vals(wrk_section,"MIXED",error=error) 00972 CALL section_vals_get(colvar%u_param%mixed_energy_section,explicit=use_mixed_energy,error=error) 00973 IF (.NOT.use_mixed_energy) NULLIFY(colvar%u_param%mixed_energy_section) 00974 ELSE IF (my_subsection(20)) THEN 00975 ! Wc hydrogen bond 00976 wrk_section => Wc_section 00977 CALL colvar_create(colvar, Wc_colvar_id, error) 00978 CALL colvar_check_points(colvar, Wc_section, error) 00979 CALL section_vals_val_get(Wc_section,"ATOMS",i_vals=iatms,error=error) 00980 CALL section_vals_val_get(wrk_section,"RCUT",r_val=my_val(1),error=error) 00981 colvar%Wc%rcut = cp_unit_to_cp2k(my_val(1),"angstrom",error=error) 00982 colvar%Wc%ids = iatms 00983 ELSE IF (my_subsection(21)) THEN 00984 ! HBP colvar 00985 wrk_section => HBP_section 00986 CALL colvar_create(colvar, HBP_colvar_id, error) 00987 CALL colvar_check_points(colvar, HBP_section, error) 00988 CALL section_vals_val_get(wrk_section,"NPOINTS",i_val=colvar%HBP%nPoints,error=error) 00989 CALL section_vals_val_get(wrk_section,"RCUT",r_val=my_val(1),error=error) 00990 colvar%HBP%rcut = cp_unit_to_cp2k(my_val(1),"angstrom",error=error) 00991 CALL section_vals_val_get(wrk_section,"RCUT",r_val=colvar%HBP%shift,error=error) 00992 00993 ALLOCATE(colvar%HBP%ids(colvar%HBP%nPoints,3),STAT=stat) 00994 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00995 ALLOCATE(colvar%HBP%ewc(colvar%HBP%nPoints),STAT=stat) 00996 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 00997 DO i = 1, colvar%HBP%nPoints 00998 CALL section_vals_val_get(wrk_section,"ATOMS",i_rep_val=i,i_vals=iatms,error=error) 00999 colvar%HBP%ids(i,:) = iatms 01000 ENDDO 01001 ELSE IF (my_subsection(22)) THEN 01002 ! Ring Puckering 01003 CALL colvar_create(colvar,ring_puckering_colvar_id, error) 01004 CALL section_vals_val_get(ring_puckering_section,"ATOMS",i_vals=iatms,error=error) 01005 colvar%ring_puckering_param%nring = SIZE(iatms) 01006 ALLOCATE(colvar%ring_puckering_param%atoms(SIZE(iatms)),stat=stat) 01007 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 01008 colvar%ring_puckering_param%atoms = iatms 01009 CALL section_vals_val_get(ring_puckering_section,"COORDINATE",& 01010 i_val=colvar%ring_puckering_param%iq,error=error) 01011 ! test the validity of the parameters 01012 ndim = colvar%ring_puckering_param%nring 01013 CALL cp_assert(ndim > 3,cp_fatal_level,cp_assertion_failed,routineP,& 01014 "CV Ring Puckering: Ring size has to be 4 or larger. "//& 01015 CPSourceFileRef) 01016 ii = colvar%ring_puckering_param%iq 01017 CALL cp_assert( (ABS(ii)/=1 .AND. ii>=-(ndim-1)/2 .AND. ii<=ndim/2) ,& 01018 cp_fatal_level,cp_assertion_failed,routineP,& 01019 "CV Ring Puckering: Invalid coordinate number. "//& 01020 CPSourceFileRef) 01021 ELSE IF (my_subsection(23)) THEN 01022 ! Minimum Distance 01023 wrk_section => mindist_section 01024 CALL colvar_create(colvar, mindist_colvar_id, error) 01025 CALL colvar_check_points(colvar, mindist_section, error) 01026 NULLIFY(colvar%mindist_param%i_dist_from,colvar%mindist_param%i_coord_from,& 01027 colvar%mindist_param%k_coord_from,colvar%mindist_param%i_coord_to,& 01028 colvar%mindist_param%k_coord_to) 01029 CALL section_vals_val_get(mindist_section,"ATOMS_DISTANCE",i_vals=iatms,error=error) 01030 colvar%mindist_param%n_dist_from = SIZE(iatms) 01031 ALLOCATE(colvar%mindist_param%i_dist_from(SIZE(iatms)),stat=stat) 01032 colvar%mindist_param%i_dist_from = iatms 01033 CALL section_vals_val_get(mindist_section,"ATOMS_FROM",n_rep_val=n_var,error=error) 01034 ndim = 0 01035 IF (n_var /= 0) THEN 01036 ! INDEX LIST 01037 DO k = 1, n_var 01038 CALL section_vals_val_get(mindist_section,"ATOMS_FROM",i_rep_val=k,i_vals=iatms,error=error) 01039 CALL reallocate(colvar%mindist_param%i_coord_from,1, ndim+SIZE(iatms)) 01040 colvar%mindist_param%i_coord_from(ndim+1:ndim+SIZE(iatms)) = iatms 01041 ndim = ndim + SIZE(iatms) 01042 END DO 01043 colvar%mindist_param%n_coord_from = ndim 01044 colvar%mindist_param%use_kinds_from = .FALSE. 01045 ELSE 01046 !KINDS 01047 CALL section_vals_val_get(mindist_section,"KINDS_FROM",n_rep_val=n_var,error=error) 01048 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 01049 DO k = 1, n_var 01050 CALL section_vals_val_get(mindist_section,"KINDS_FROM",i_rep_val=k,c_vals=c_kinds,error=error) 01051 CALL reallocate(colvar%mindist_param%k_coord_from,1, ndim+SIZE(c_kinds)) 01052 colvar%mindist_param%k_coord_from(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 01053 ndim = ndim + SIZE(c_kinds) 01054 END DO 01055 colvar%mindist_param%n_coord_from = 0 01056 colvar%mindist_param%use_kinds_from = .TRUE. 01057 ! Uppercase the label 01058 DO k = 1, ndim 01059 CALL uppercase(colvar%mindist_param%k_coord_from(k)) 01060 END DO 01061 END IF 01062 01063 01064 CALL section_vals_val_get(mindist_section,"ATOMS_TO",n_rep_val=n_var,error=error) 01065 ndim = 0 01066 IF (n_var /= 0) THEN 01067 ! INDEX LIST 01068 DO k = 1, n_var 01069 CALL section_vals_val_get(mindist_section,"ATOMS_TO",i_rep_val=k,i_vals=iatms,error=error) 01070 CALL reallocate(colvar%mindist_param%i_coord_to,1, ndim+SIZE(iatms)) 01071 colvar%mindist_param%i_coord_to(ndim+1:ndim+SIZE(iatms)) = iatms 01072 ndim = ndim + SIZE(iatms) 01073 END DO 01074 colvar%mindist_param%n_coord_to = ndim 01075 colvar%mindist_param%use_kinds_to = .FALSE. 01076 ELSE 01077 !KINDS 01078 CALL section_vals_val_get(mindist_section,"KINDS_TO",n_rep_val=n_var,error=error) 01079 CPPostcondition(n_var>0,cp_failure_level,routinep,error,failure) 01080 DO k = 1, n_var 01081 CALL section_vals_val_get(mindist_section,"KINDS_TO",i_rep_val=k,c_vals=c_kinds,error=error) 01082 CALL reallocate(colvar%mindist_param%k_coord_to,1, ndim+SIZE(c_kinds)) 01083 colvar%mindist_param%k_coord_to(ndim+1:ndim+SIZE(c_kinds)) = c_kinds 01084 ndim = ndim + SIZE(c_kinds) 01085 END DO 01086 colvar%mindist_param%n_coord_to = 0 01087 colvar%mindist_param%use_kinds_to = .TRUE. 01088 ! Uppercase the label 01089 DO k = 1, ndim 01090 CALL uppercase(colvar%mindist_param%k_coord_to(k)) 01091 END DO 01092 END IF 01093 01094 CALL section_vals_val_get(mindist_section,"R0",r_val=colvar%mindist_param%r_cut,error=error) 01095 CALL section_vals_val_get(mindist_section,"NN",i_val=colvar%mindist_param%p_exp,error=error) 01096 CALL section_vals_val_get(mindist_section,"ND",i_val=colvar%mindist_param%q_exp,error=error) 01097 ! CALL section_vals_val_get(mindist_section,"NC",r_val=colvar%mindist_param%n_cut,error=error) 01098 CALL section_vals_val_get(mindist_section,"LAMBDA",r_val=colvar%mindist_param%lambda,error=error) 01099 END IF 01100 CALL colvar_setup(colvar, error) 01101 01102 iw = cp_print_key_unit_nr(logger,colvar_section,& 01103 "PRINT%PROGRAM_RUN_INFO",extension=".colvarLog",error=error) 01104 IF (iw>0) THEN 01105 tag = "ATOMS: " 01106 IF (colvar%use_points) tag = "POINTS:" 01107 ! Description header 01108 IF (colvar%type_id/=combine_colvar_id) THEN 01109 WRITE ( iw, '( A )')' '//& 01110 '----------------------------------------------------------------------' 01111 WRITE ( iw, '( A,I8)' ) ' COLVARS| COLVAR INPUT INDEX: ',icol 01112 END IF 01113 ! Colvar Description 01114 SELECT CASE(colvar%type_id) 01115 CASE(angle_colvar_id) 01116 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| ANGLE >>> '//tag,& 01117 colvar%angle_param%i_at_angle 01118 CASE(dfunct_colvar_id) 01119 WRITE ( iw, '( A,T49,4I8)' ) ' COLVARS| DISTANCE DIFFERENCE >>> '//tag,& 01120 colvar%dfunct_param%i_at_dfunct 01121 CASE(plane_distance_colvar_id) 01122 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| PLANE DISTANCE - PLANE >>> '//tag,& 01123 colvar%plane_distance_param%plane 01124 WRITE ( iw, '( A,T73,1I8)' ) ' COLVARS| PLANE DISTANCE - POINT >>> '//tag,& 01125 colvar%plane_distance_param%point 01126 CASE(plane_plane_angle_colvar_id) 01127 IF (colvar%plane_plane_angle_param%plane1%type_of_def==plane_def_atoms) THEN 01128 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag,& 01129 colvar%plane_plane_angle_param%plane1%points 01130 ELSE 01131 WRITE ( iw, '( A,T57,3F8.3)' ) ' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag,& 01132 colvar%plane_plane_angle_param%plane1%normal_vec 01133 END IF 01134 01135 IF (colvar%plane_plane_angle_param%plane2%type_of_def==plane_def_atoms) THEN 01136 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (ATOMS) >>> '//tag,& 01137 colvar%plane_plane_angle_param%plane2%points 01138 ELSE 01139 WRITE ( iw, '( A,T57,3F8.3)' ) ' COLVARS| PLANE-PLANE ANGLE - PLANE 1 (VECTOR) >>> '//tag,& 01140 colvar%plane_plane_angle_param%plane2%normal_vec 01141 END IF 01142 CASE(torsion_colvar_id) 01143 WRITE ( iw, '( A,T49,4I8)' ) ' COLVARS| TORSION >>> '//tag,& 01144 colvar%torsion_param%i_at_tors 01145 CASE(dist_colvar_id) 01146 WRITE ( iw, '( A,T65,2I8)' ) ' COLVARS| BOND >>> '//tag,& 01147 colvar%dist_param%i_at,colvar%dist_param%j_at 01148 CASE(coord_colvar_id) 01149 IF(colvar%coord_param%do_chain) THEN 01150 WRITE ( iw, '( A)' ) ' COLVARS| COORDINATION CHAIN FC(from->to)*FC(to->to_B)>> ' 01151 END IF 01152 IF (colvar%coord_param%use_kinds_from) THEN 01153 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| COORDINATION >>> FROM KINDS',& 01154 ADJUSTR(colvar%coord_param%c_kinds_from(kk)(1:10)),& 01155 kk=1,SIZE(colvar%coord_param%c_kinds_from)) 01156 ELSE 01157 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COORDINATION >>> FROM '//tag,& 01158 colvar%coord_param%i_at_from(kk),& 01159 kk=1,SIZE(colvar%coord_param%i_at_from)) 01160 END IF 01161 IF (colvar%coord_param%use_kinds_to) THEN 01162 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| COORDINATION >>> TO KINDS',& 01163 ADJUSTR(colvar%coord_param%c_kinds_to(kk)(1:10)),& 01164 kk=1,SIZE(colvar%coord_param%c_kinds_to)) 01165 ELSE 01166 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COORDINATION >>> TO '//tag,& 01167 colvar%coord_param%i_at_to(kk),& 01168 kk=1,SIZE(colvar%coord_param%i_at_to)) 01169 END IF 01170 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| R0',colvar%coord_param%r_0 01171 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| NN',colvar%coord_param%nncrd 01172 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| ND',colvar%coord_param%ndcrd 01173 IF(colvar%coord_param%do_chain) THEN 01174 IF (colvar%coord_param%use_kinds_to_b) THEN 01175 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| COORDINATION >>> TO KINDS B',& 01176 ADJUSTR(colvar%coord_param%c_kinds_to_b(kk)(1:10)),& 01177 kk=1,SIZE(colvar%coord_param%c_kinds_to_b)) 01178 ELSE 01179 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COORDINATION >>> TO '//tag//' B',& 01180 colvar%coord_param%i_at_to_b(kk),& 01181 kk=1,SIZE(colvar%coord_param%i_at_to_b)) 01182 END IF 01183 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| R0 B',colvar%coord_param%r_0_b 01184 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| NN B',colvar%coord_param%nncrd_b 01185 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| ND B',colvar%coord_param%ndcrd_b 01186 END IF 01187 CASE(population_colvar_id) 01188 IF (colvar%population_param%use_kinds_from) THEN 01189 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| POPULATION based on coordination >>> FROM KINDS',& 01190 ADJUSTR(colvar%population_param%c_kinds_from(kk)(1:10)),& 01191 kk=1,SIZE(colvar%population_param%c_kinds_from)) 01192 ELSE 01193 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| POPULATION based on coordination >>> FROM '//tag,& 01194 colvar%population_param%i_at_from(kk),& 01195 kk=1,SIZE(colvar%population_param%i_at_from)) 01196 END IF 01197 IF (colvar%population_param%use_kinds_to) THEN 01198 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| POPULATION based on coordination >>> TO KINDS',& 01199 ADJUSTR(colvar%population_param%c_kinds_to(kk)(1:10)),& 01200 kk=1,SIZE(colvar%population_param%c_kinds_to)) 01201 ELSE 01202 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| POPULATION based on coordination >>> TO '//tag,& 01203 colvar%population_param%i_at_to(kk),& 01204 kk=1,SIZE(colvar%population_param%i_at_to)) 01205 END IF 01206 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| R0',colvar%population_param%r_0 01207 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| NN',colvar%population_param%nncrd 01208 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| ND',colvar%population_param%ndcrd 01209 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| N0',colvar%population_param%n0 01210 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| SIGMA',colvar%population_param%sigma 01211 CASE(gyration_colvar_id) 01212 IF (colvar%gyration_param%use_kinds) THEN 01213 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| Gyration Radius >>> KINDS',& 01214 ADJUSTR(colvar%gyration_param%c_kinds(kk)(1:10)),& 01215 kk=1,SIZE(colvar%gyration_param%c_kinds)) 01216 ELSE 01217 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| Gyration Radius >>> ATOMS '//tag,& 01218 colvar%gyration_param%i_at(kk),& 01219 kk=1,SIZE(colvar%gyration_param%i_at)) 01220 END IF 01221 CASE(rotation_colvar_id) 01222 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| BOND_ROTATION - POINT 1 LINE 1 >>> '//tag,& 01223 colvar%rotation_param%i_at1_bond1 01224 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| BOND_ROTATION - POINT 2 LINE 1 >>> '//tag,& 01225 colvar%rotation_param%i_at2_bond1 01226 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| BOND_ROTATION - POINT 1 LINE 2 >>> '//tag,& 01227 colvar%rotation_param%i_at1_bond2 01228 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| BOND_ROTATION - POINT 2 LINE 2 >>> '//tag,& 01229 colvar%rotation_param%i_at2_bond2 01230 CASE(qparm_colvar_id) 01231 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| Q-PARM >>> FROM '//tag,& 01232 colvar%qparm_param%i_at_from(kk),& 01233 kk=1,SIZE(colvar%qparm_param%i_at_from)) 01234 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| Q-PARM >>> TO '//tag,& 01235 colvar%qparm_param%i_at_to(kk),& 01236 kk=1,SIZE(colvar%qparm_param%i_at_to)) 01237 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| RCUT',colvar%qparm_param%rcut 01238 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| ALPHA',colvar%qparm_param%alpha 01239 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| L',colvar%qparm_param%l 01240 CASE(combine_colvar_id) 01241 WRITE ( iw, '( A)' ) ' COLVARS| COMBINING FUNCTION : '//& 01242 TRIM(colvar%combine_cvs_param%function) 01243 WRITE ( iw, '( A)', ADVANCE="NO") ' COLVARS| VARIABLES : ' 01244 DO i = 1, SIZE(colvar%combine_cvs_param%variables) 01245 WRITE ( iw, '( A)', ADVANCE="NO")& 01246 TRIM(colvar%combine_cvs_param%variables(i))//" " 01247 END DO 01248 WRITE ( iw, '(/)') 01249 WRITE ( iw, '( A)' ) ' COLVARS| DEFINED PARAMETERS [label] [value]:' 01250 DO i = 1, SIZE(colvar%combine_cvs_param%c_parameters) 01251 WRITE ( iw, '( A,A7,F9.3)' ) ' ',& 01252 TRIM(colvar%combine_cvs_param%c_parameters(i)),colvar%combine_cvs_param%v_parameters(i) 01253 END DO 01254 WRITE ( iw, '( A,T71,G10.5)' ) ' COLVARS| ERROR ON DERIVATIVE EVALUATION',& 01255 colvar%combine_cvs_param%lerr 01256 WRITE ( iw, '( A,T71,G10.5)' ) ' COLVARS| DX',& 01257 colvar%combine_cvs_param%dx 01258 CASE(reaction_path_colvar_id) 01259 CALL cp_unimplemented_error(fromWhere=routineP, & 01260 message="Description header for REACTION_PATH COLVAR missing!!", & 01261 error=error, error_level=cp_warning_level) 01262 CASE(distance_from_path_colvar_id) 01263 CALL cp_unimplemented_error(fromWhere=routineP, & 01264 message="Description header for REACTION_PATH COLVAR missing!!", & 01265 error=error, error_level=cp_warning_level) 01266 CASE(hydronium_colvar_id) 01267 CALL cp_unimplemented_error(fromWhere=routineP, & 01268 message="Description header for HYDRONIUM COLVAR missing!!", & 01269 error=error, error_level=cp_warning_level) 01270 CASE(rmsd_colvar_id) 01271 CALL cp_unimplemented_error(fromWhere=routineP, & 01272 message="Description header for RMSD COLVAR missing!!", & 01273 error=error, error_level=cp_warning_level) 01274 CASE(xyz_diag_colvar_id) 01275 NULLIFY(section, keyword, enum) 01276 CALL create_colvar_xyz_d_section(section,error=error) 01277 keyword => section_get_keyword(section,"COMPONENT",error=error) 01278 CALL keyword_get(keyword,enum=enum,error=error) 01279 tag_comp = enum_i2c(enum,colvar%xyz_diag_param%component,error=error) 01280 CALL section_release(section,error=error) 01281 01282 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| POSITION ('//TRIM(tag_comp)& 01283 //') >>> '//tag,colvar%xyz_diag_param%i_atom 01284 CASE(xyz_outerdiag_colvar_id) 01285 NULLIFY(section, keyword, enum) 01286 CALL create_colvar_xyz_od_section(section,error=error) 01287 keyword => section_get_keyword(section,"COMPONENT_A",error=error) 01288 CALL keyword_get(keyword,enum=enum,error=error) 01289 tag_comp1 = enum_i2c(enum,colvar%xyz_outerdiag_param%components(1),error=error) 01290 keyword => section_get_keyword(section,"COMPONENT_B",error=error) 01291 CALL keyword_get(keyword,enum=enum,error=error) 01292 tag_comp2 = enum_i2c(enum,colvar%xyz_outerdiag_param%components(2),error=error) 01293 CALL section_release(section,error=error) 01294 01295 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| CROSS TERM POSITION ('//TRIM(tag_comp1)& 01296 //" * "//TRIM(tag_comp2)//') >>> '//tag, colvar%xyz_outerdiag_param%i_atoms 01297 CASE(u_colvar_id) 01298 WRITE ( iw, '( A,T77,A4)' ) ' COLVARS| ENERGY >>> '//tag,'all!' 01299 CASE(Wc_colvar_id) 01300 WRITE ( iw, '( A,T57,F16.8)' ) ' COLVARS| Wc >>> RCUT: ',& 01301 colvar%Wc%rcut 01302 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| Wc >>> '//tag,& 01303 colvar%Wc%ids 01304 CASE(HBP_colvar_id) 01305 WRITE ( iw, '( A,T57,I8)' ) ' COLVARS| HBP >>> NPOINTS',& 01306 colvar%HBP%nPoints 01307 WRITE ( iw, '( A,T57,F16.8)' ) ' COLVARS| HBP >>> RCUT',& 01308 colvar%HBP%rcut 01309 WRITE ( iw, '( A,T57,F16.8)' ) ' COLVARS| HBP >>> RCUT',& 01310 colvar%HBP%shift 01311 DO i =1, colvar%HBP%nPoints 01312 WRITE ( iw, '( A,T57,3I8)' ) ' COLVARS| HBP >>> '//tag,& 01313 colvar%HBP%ids(i,:) 01314 ENDDO 01315 CASE(ring_puckering_colvar_id) 01316 WRITE ( iw, '( A,T57,I8)' ) ' COLVARS| Ring Puckering >>> ring size',& 01317 colvar%ring_puckering_param%nring 01318 IF (colvar%ring_puckering_param%iq == 0) THEN 01319 WRITE ( iw, '( A,T40,A)' ) ' COLVARS| Ring Puckering >>> coordinate',& 01320 ' Total Puckering Amplitude' 01321 ELSEIF (colvar%ring_puckering_param%iq > 0) THEN 01322 WRITE ( iw, '( A,T35,A,T57,I8)' ) ' COLVARS| Ring Puckering >>> coordinate',& 01323 ' Puckering Amplitude',& 01324 colvar%ring_puckering_param%iq 01325 ELSE 01326 WRITE ( iw, '( A,T35,A,T57,I8)' ) ' COLVARS| Ring Puckering >>> coordinate',& 01327 ' Puckering Angle',& 01328 colvar%ring_puckering_param%iq 01329 END IF 01330 CASE(mindist_colvar_id) 01331 WRITE ( iw, '( A)' ) ' COLVARS| CONDITIONED DISTANCE>> ' 01332 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COND.DISTANCE >>> DISTANCE FROM '//tag,& 01333 colvar%mindist_param%i_dist_from(kk),& 01334 kk=1,SIZE(colvar%mindist_param%i_dist_from)) 01335 IF (colvar%mindist_param%use_kinds_from) THEN 01336 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| COND.DIST. >>> COORDINATION FROM KINDS ',& 01337 ADJUSTR(colvar%mindist_param%k_coord_from(kk)(1:10)),& 01338 kk=1,SIZE(colvar%mindist_param%k_coord_from)) 01339 ELSE 01340 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COND.DIST. >>> COORDINATION FROM '//tag,& 01341 colvar%mindist_param%i_coord_from(kk),& 01342 kk=1,SIZE(colvar%mindist_param%i_coord_from)) 01343 END IF 01344 IF (colvar%mindist_param%use_kinds_to) THEN 01345 WRITE ( iw, '( A,T71,A10)' ) (' COLVARS| COND.DIST. >>> COORDINATION TO KINDS ',& 01346 ADJUSTR(colvar%mindist_param%k_coord_to(kk)(1:10)),& 01347 kk=1,SIZE(colvar%mindist_param%k_coord_to)) 01348 ELSE 01349 WRITE ( iw, '( A,T71,I10)' ) (' COLVARS| COND.DIST. >>> COORDINATION TO '//tag,& 01350 colvar%mindist_param%i_coord_to(kk),& 01351 kk=1,SIZE(colvar%mindist_param%i_coord_to)) 01352 END IF 01353 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| R0',colvar%mindist_param%r_cut 01354 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| NN',colvar%mindist_param%p_exp 01355 WRITE ( iw, '( A,T71,I10)' ) ' COLVARS| ND',colvar%mindist_param%q_exp 01356 WRITE ( iw, '( A,T71,F10.5)' ) ' COLVARS| LAMBDA',colvar%mindist_param%lambda 01357 01358 END SELECT 01359 IF (colvar%use_points) THEN 01360 WRITE ( iw, '( A)') ' COLVARS| INFORMATION ON DEFINED GEOMETRICAL POINTS' 01361 DO kk = 1, SIZE(colvar%points) 01362 point_section => section_vals_get_subs_vals(wrk_section,"POINT",error=error) 01363 CALL section_vals_val_get(point_section,"TYPE",i_rep_section=kk,c_val=tmpStr,error=error) 01364 tmpStr2 = cp_to_string(kk) 01365 WRITE ( iw, '( A)') ' COLVARS| POINT Nr.'//TRIM(tmpStr2)//' OF TYPE: '//TRIM(tmpStr) 01366 IF (ASSOCIATED(colvar%points(kk)%atoms)) THEN 01367 WRITE ( iw, '( A)') ' COLVARS| ATOMS BUILDING THE GEOMETRICAL POINT' 01368 WRITE ( iw, '( A, I10)')(' COLVARS| ATOM:',colvar%points(kk)%atoms(k),k=1,SIZE(colvar%points(kk)%atoms)) 01369 ELSE 01370 WRITE ( iw, '( A,4X,3F12.6)') ' COLVARS| XYZ POSITION OF FIXED POINT:', colvar%points(kk)%r 01371 END IF 01372 END DO 01373 END IF 01374 ! Close the description layer 01375 IF (colvar%type_id/=combine_colvar_id) THEN 01376 WRITE ( iw, '( A )')' '//& 01377 '----------------------------------------------------------------------' 01378 ELSE 01379 WRITE ( iw, '( A )')' '//& 01380 '**********************************************************************' 01381 END IF 01382 END IF 01383 CALL cp_print_key_finished_output(iw,logger,colvar_section,& 01384 "PRINT%PROGRAM_RUN_INFO", error=error) 01385 CALL timestop(handle) 01386 END SUBROUTINE colvar_read 01387 01388 ! ***************************************************************************** 01396 SUBROUTINE colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list, error) 01397 TYPE(colvar_type), POINTER :: colvar 01398 TYPE(cell_type), POINTER :: cell 01399 TYPE(particle_type), DIMENSION(:), 01400 POINTER :: particles 01401 REAL(kind=dp), DIMENSION(:, :), 01402 INTENT(IN), OPTIONAL :: pos 01403 TYPE(fixd_constraint_type), 01404 DIMENSION(:), OPTIONAL, POINTER :: fixd_list 01405 TYPE(cp_error_type), INTENT(inout) :: error 01406 01407 CHARACTER(len=*), PARAMETER :: routineN = 'colvar_eval_mol_f', 01408 routineP = moduleN//':'//routineN 01409 01410 INTEGER :: i, j 01411 LOGICAL :: colvar_ok, failure 01412 01413 failure=.FALSE. 01414 01415 colvar_ok=ASSOCIATED(colvar) 01416 CPAssert(colvar_ok,cp_failure_level,routineP,error,failure) 01417 IF (.NOT. failure) THEN 01418 IF (PRESENT(pos)) THEN 01419 DO i = 1, SIZE(colvar%i_atom) 01420 j = colvar%i_atom(i) 01421 particles(j)%r=pos(:,j) 01422 END DO 01423 END IF 01424 ! Initialize the content of the derivative 01425 colvar%dsdr=0.0_dp 01426 SELECT CASE(colvar%type_id) 01427 CASE (dist_colvar_id) 01428 CALL dist_colvar(colvar,cell,particles=particles,error=error) 01429 CASE (coord_colvar_id) 01430 CALL coord_colvar(colvar,cell,particles=particles, error=error) 01431 CASE (population_colvar_id) 01432 CALL population_colvar(colvar,cell,particles=particles, error=error) 01433 CASE (gyration_colvar_id) 01434 CALL gyration_radius_colvar(colvar,cell,particles=particles, error=error) 01435 CASE (torsion_colvar_id) 01436 CALL torsion_colvar(colvar, cell, particles=particles, error=error) 01437 CASE (angle_colvar_id) 01438 CALL angle_colvar(colvar, cell, particles=particles, error=error) 01439 CASE (dfunct_colvar_id) 01440 CALL dfunct_colvar(colvar, cell, particles=particles, error=error) 01441 CASE (plane_distance_colvar_id) 01442 CALL plane_distance_colvar(colvar,cell, particles=particles, error=error) 01443 CASE (plane_plane_angle_colvar_id) 01444 CALL plane_plane_angle_colvar(colvar,cell, particles=particles, error=error) 01445 CASE (rotation_colvar_id) 01446 CALL rotation_colvar(colvar,cell, particles=particles, error=error) 01447 CASE (qparm_colvar_id) 01448 CALL qparm_colvar(colvar,cell, particles=particles, error=error) 01449 CASE (hydronium_colvar_id) 01450 CALL hydronium_colvar(colvar,cell,particles=particles, error=error) 01451 CASE(rmsd_colvar_id) 01452 CALL rmsd_colvar(colvar,particles=particles,error=error) 01453 CASE (reaction_path_colvar_id) 01454 CALL reaction_path_colvar(colvar,cell,particles=particles, error=error) 01455 CASE (distance_from_path_colvar_id) 01456 CALL distance_from_path_colvar(colvar,cell,particles=particles, error=error) 01457 CASE (combine_colvar_id) 01458 CALL combine_colvar(colvar,cell,particles=particles, error=error) 01459 CASE (xyz_diag_colvar_id) 01460 CALL xyz_diag_colvar(colvar,cell,particles=particles, error=error) 01461 CASE (xyz_outerdiag_colvar_id) 01462 CALL xyz_outerdiag_colvar(colvar,cell,particles=particles, error=error) 01463 CASE (ring_puckering_colvar_id) 01464 CALL ring_puckering_colvar(colvar,cell,particles=particles, error=error) 01465 CASE (mindist_colvar_id) 01466 CALL mindist_colvar(colvar,cell,particles=particles, error=error) 01467 CASE (u_colvar_id) 01468 CALL cp_unimplemented_error(fromWhere=routineP, & 01469 message="need force_env!", error=error, error_level=cp_failure_level) 01470 CASE (Wc_colvar_id) 01471 !!! FIXME this is rubbish at the moment as we have no force to be computed on this 01472 CALL Wc_colvar(colvar, cell, particles=particles, error=error) 01473 CASE (HBP_colvar_id) 01474 !!! FIXME this is rubbish at the moment as we have no force to be computed on this 01475 CALL HBP_colvar(colvar, cell, particles=particles, error=error) 01476 CASE DEFAULT 01477 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01478 END SELECT 01479 ! Check for fixed atom constraints 01480 IF (PRESENT(fixd_list)) CALL check_fixed_atom_cns_colv(fixd_list, colvar) 01481 END IF 01482 01483 END SUBROUTINE colvar_eval_mol_f 01484 01485 ! ***************************************************************************** 01495 SUBROUTINE colvar_eval_glob_f(icolvar,force_env,error) 01496 INTEGER :: icolvar 01497 TYPE(force_env_type), POINTER :: force_env 01498 TYPE(cp_error_type), INTENT(inout) :: error 01499 01500 CHARACTER(len=*), PARAMETER :: routineN = 'colvar_eval_glob_f', 01501 routineP = moduleN//':'//routineN 01502 01503 LOGICAL :: colvar_ok, failure 01504 TYPE(cell_type), POINTER :: cell 01505 TYPE(colvar_type), POINTER :: colvar 01506 TYPE(cp_subsys_type), POINTER :: subsys 01507 TYPE(qs_environment_type), POINTER :: qs_env 01508 01509 failure=.FALSE. 01510 NULLIFY(subsys,cell,colvar,qs_env) 01511 CALL force_env_get(force_env,subsys=subsys,cell=cell,qs_env=qs_env,error=error) 01512 colvar_ok=ASSOCIATED(subsys%colvar_p) 01513 CPAssert(colvar_ok,cp_failure_level,routineP,error,failure) 01514 01515 IF (.NOT. failure) THEN 01516 colvar => subsys%colvar_p(icolvar)%colvar 01517 ! Initialize the content of the derivative 01518 colvar%dsdr=0.0_dp 01519 SELECT CASE(colvar%type_id) 01520 CASE (dist_colvar_id) 01521 CALL dist_colvar(colvar,cell,subsys=subsys,error=error) 01522 CASE (coord_colvar_id) 01523 CALL coord_colvar(colvar,cell,subsys=subsys, error=error) 01524 CASE (population_colvar_id) 01525 CALL population_colvar(colvar,cell,subsys=subsys, error=error) 01526 CASE (gyration_colvar_id) 01527 CALL gyration_radius_colvar(colvar,cell,subsys=subsys, error=error) 01528 CASE (torsion_colvar_id) 01529 CALL torsion_colvar(colvar,cell,subsys=subsys, no_riemann_sheet_op=.TRUE.,error=error) 01530 CASE (angle_colvar_id) 01531 CALL angle_colvar(colvar,cell,subsys=subsys, error=error) 01532 CASE (dfunct_colvar_id) 01533 CALL dfunct_colvar(colvar,cell,subsys=subsys, error=error) 01534 CASE (plane_distance_colvar_id) 01535 CALL plane_distance_colvar(colvar,cell,subsys=subsys, error=error) 01536 CASE (plane_plane_angle_colvar_id) 01537 CALL plane_plane_angle_colvar(colvar,cell,subsys=subsys, error=error) 01538 CASE (rotation_colvar_id) 01539 CALL rotation_colvar(colvar,cell,subsys=subsys, error=error) 01540 CASE (qparm_colvar_id) 01541 CALL qparm_colvar(colvar,cell,subsys=subsys, error=error) 01542 CASE (hydronium_colvar_id) 01543 CALL hydronium_colvar(colvar,cell,subsys=subsys, error=error) 01544 CASE(rmsd_colvar_id) 01545 CALL rmsd_colvar(colvar,subsys=subsys,error=error) 01546 CASE (reaction_path_colvar_id) 01547 CALL reaction_path_colvar(colvar,cell,subsys=subsys, error=error) 01548 CASE (distance_from_path_colvar_id) 01549 CALL distance_from_path_colvar(colvar,cell, subsys=subsys, error=error) 01550 CASE (combine_colvar_id) 01551 CALL combine_colvar(colvar,cell,subsys=subsys, error=error) 01552 CASE (xyz_diag_colvar_id) 01553 CALL xyz_diag_colvar(colvar,cell,subsys=subsys, error=error) 01554 CASE (xyz_outerdiag_colvar_id) 01555 CALL xyz_outerdiag_colvar(colvar,cell,subsys=subsys, error=error) 01556 CASE (u_colvar_id) 01557 CALL u_colvar(colvar,force_env=force_env,subsys=subsys,error=error) 01558 CASE (Wc_colvar_id) 01559 CALL Wc_colvar(colvar,cell,subsys=subsys, qs_env=qs_env, error=error) 01560 CASE (HBP_colvar_id) 01561 CALL HBP_colvar(colvar,cell,subsys=subsys, qs_env=qs_env, error=error) 01562 CASE (ring_puckering_colvar_id) 01563 CALL ring_puckering_colvar(colvar,cell,subsys=subsys, error=error) 01564 CASE (mindist_colvar_id) 01565 CALL mindist_colvar(colvar,cell,subsys=subsys , error=error) 01566 CASE DEFAULT 01567 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01568 END SELECT 01569 ! Check for fixed atom constraints 01570 CALL check_fixed_atom_cns_colv(subsys%gci%fixd_list,colvar) 01571 END IF 01572 END SUBROUTINE colvar_eval_glob_f 01573 01574 ! ***************************************************************************** 01582 SUBROUTINE colvar_recursive_eval(colvar,cell,particles,error) 01583 TYPE(colvar_type), POINTER :: colvar 01584 TYPE(cell_type), POINTER :: cell 01585 TYPE(particle_type), DIMENSION(:), 01586 POINTER :: particles 01587 TYPE(cp_error_type), INTENT(inout) :: error 01588 01589 CHARACTER(len=*), PARAMETER :: routineN = 'colvar_recursive_eval', 01590 routineP = moduleN//':'//routineN 01591 01592 LOGICAL :: failure 01593 01594 failure=.FALSE. 01595 01596 IF (.NOT. failure) THEN 01597 ! Initialize the content of the derivative 01598 colvar%dsdr=0.0_dp 01599 SELECT CASE(colvar%type_id) 01600 CASE (dist_colvar_id) 01601 CALL dist_colvar(colvar,cell,particles=particles,error=error) 01602 CASE (coord_colvar_id) 01603 CALL coord_colvar(colvar,cell,particles=particles, error=error) 01604 CASE (torsion_colvar_id) 01605 CALL torsion_colvar(colvar, cell, particles=particles, error=error) 01606 CASE (angle_colvar_id) 01607 CALL angle_colvar(colvar, cell, particles=particles, error=error) 01608 CASE (dfunct_colvar_id) 01609 CALL dfunct_colvar(colvar, cell, particles=particles, error=error) 01610 CASE (plane_distance_colvar_id) 01611 CALL plane_distance_colvar(colvar,cell, particles=particles, error=error) 01612 CASE (plane_plane_angle_colvar_id) 01613 CALL plane_plane_angle_colvar(colvar,cell, particles=particles, error=error) 01614 CASE (rotation_colvar_id) 01615 CALL rotation_colvar(colvar,cell, particles=particles, error=error) 01616 CASE (qparm_colvar_id) 01617 CALL qparm_colvar(colvar,cell, particles=particles, error=error) 01618 CASE (hydronium_colvar_id) 01619 CALL hydronium_colvar(colvar,cell,particles=particles, error=error) 01620 CASE(rmsd_colvar_id) 01621 CALL rmsd_colvar(colvar,particles=particles,error=error) 01622 CASE (reaction_path_colvar_id) 01623 CALL reaction_path_colvar(colvar,cell,particles=particles, error=error) 01624 CASE (distance_from_path_colvar_id) 01625 CALL distance_from_path_colvar(colvar,cell,particles=particles, error=error) 01626 CASE (combine_colvar_id) 01627 CALL combine_colvar(colvar,cell,particles=particles, error=error) 01628 CASE (xyz_diag_colvar_id) 01629 CALL xyz_diag_colvar(colvar,cell,particles=particles, error=error) 01630 CASE (xyz_outerdiag_colvar_id) 01631 CALL xyz_outerdiag_colvar(colvar,cell,particles=particles, error=error) 01632 CASE (ring_puckering_colvar_id) 01633 CALL ring_puckering_colvar(colvar,cell,particles=particles, error=error) 01634 CASE (mindist_colvar_id) 01635 CALL mindist_colvar(colvar,cell,particles=particles, error=error) 01636 CASE (u_colvar_id) 01637 CALL cp_unimplemented_error(fromWhere=routineP, & 01638 message="need force_env!", error=error, error_level=cp_failure_level) 01639 CASE (Wc_colvar_id) 01640 CALL Wc_colvar(colvar, cell, particles=particles, error=error) 01641 CASE (HBP_colvar_id) 01642 CALL HBP_colvar(colvar, cell, particles=particles, error=error) 01643 CASE DEFAULT 01644 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01645 END SELECT 01646 END IF 01647 END SUBROUTINE colvar_recursive_eval 01648 01649 ! ***************************************************************************** 01653 SUBROUTINE get_coordinates(colvar, i, ri, my_particles) 01654 TYPE(colvar_type), POINTER :: colvar 01655 INTEGER, INTENT(IN) :: i 01656 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: ri 01657 TYPE(particle_type), DIMENSION(:), 01658 POINTER :: my_particles 01659 01660 IF (colvar%use_points) THEN 01661 CALL eval_point_pos(colvar%points(i),my_particles,ri) 01662 ELSE 01663 ri(:) = my_particles(i)%r(:) 01664 END IF 01665 01666 END SUBROUTINE get_coordinates 01667 01668 01669 ! ***************************************************************************** 01673 SUBROUTINE get_mass(colvar, i, mi, my_particles) 01674 TYPE(colvar_type), POINTER :: colvar 01675 INTEGER, INTENT(IN) :: i 01676 REAL(KIND=dp), INTENT(OUT) :: mi 01677 TYPE(particle_type), DIMENSION(:), 01678 POINTER :: my_particles 01679 01680 IF (colvar%use_points) THEN 01681 CALL eval_point_mass(colvar%points(i),my_particles,mi) 01682 ELSE 01683 mi = my_particles(i)%atomic_kind%mass 01684 END IF 01685 01686 END SUBROUTINE get_mass 01687 01688 ! ***************************************************************************** 01692 SUBROUTINE put_derivative(colvar, i, fi) 01693 TYPE(colvar_type), POINTER :: colvar 01694 INTEGER, INTENT(IN) :: i 01695 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: fi 01696 01697 IF (colvar%use_points) THEN 01698 CALL eval_point_der(colvar%points,i,colvar%dsdr,fi) 01699 ELSE 01700 colvar%dsdr(:,i) = colvar%dsdr(:,i) + fi 01701 END IF 01702 01703 END SUBROUTINE put_derivative 01704 01705 ! ***************************************************************************** 01709 SUBROUTINE xyz_diag_colvar(colvar,cell,subsys,particles,error) 01710 TYPE(colvar_type), POINTER :: colvar 01711 TYPE(cell_type), POINTER :: cell 01712 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 01713 TYPE(particle_type), DIMENSION(:), 01714 OPTIONAL, POINTER :: particles 01715 TYPE(cp_error_type), INTENT(inout) :: error 01716 01717 CHARACTER(len=*), PARAMETER :: routineN = 'xyz_diag_colvar', 01718 routineP = moduleN//':'//routineN 01719 01720 INTEGER :: i 01721 LOGICAL :: failure 01722 REAL(dp) :: fi(3), r, r0(3), ss(3), 01723 xi(3), xpi(3) 01724 TYPE(particle_list_type), POINTER :: particles_i 01725 TYPE(particle_type), DIMENSION(:), 01726 POINTER :: my_particles 01727 01728 failure=.FALSE. 01729 NULLIFY(particles_i) 01730 01731 CPPrecondition(colvar%type_id==xyz_diag_colvar_id,cp_failure_level,routineP,error,failure) 01732 IF (PRESENT(particles)) THEN 01733 my_particles => particles 01734 ELSE 01735 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 01736 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 01737 my_particles => particles_i%els 01738 END IF 01739 i=colvar%xyz_diag_param%i_atom 01740 ! Atom coordinates 01741 CALL get_coordinates(colvar, i, xpi, my_particles) 01742 ! Use the current coordinates as initial coordinates, if no initialization 01743 ! was performed yet 01744 IF (ALL(colvar%xyz_diag_param%r0 == HUGE(0.0_dp))) THEN 01745 colvar%xyz_diag_param%r0 = xpi 01746 END IF 01747 r0 = colvar%xyz_diag_param%r0 01748 01749 IF(colvar%xyz_diag_param%use_pbc)THEN 01750 ss=MATMUL(cell%h_inv,xpi-r0) 01751 ss=ss-NINT(ss) 01752 xi=MATMUL(cell%hmat,ss) 01753 ELSE 01754 xi=xpi-r0 01755 END IF 01756 01757 SELECT CASE (colvar%xyz_diag_param%component) 01758 CASE(do_clv_x) 01759 xi(2)=0.0_dp 01760 xi(3)=0.0_dp 01761 CASE(do_clv_y) 01762 xi(1)=0.0_dp 01763 xi(3)=0.0_dp 01764 CASE(do_clv_z) 01765 xi(1)=0.0_dp 01766 xi(2)=0.0_dp 01767 CASE(do_clv_xy) 01768 xi(3)=0.0_dp 01769 CASE(do_clv_xz) 01770 xi(2)=0.0_dp 01771 CASE(do_clv_yz) 01772 xi(1)=0.0_dp 01773 CASE DEFAULT 01774 ! do_clv_xyz 01775 END SELECT 01776 01777 r=xi(1)**2+xi(2)**2+xi(3)**2 01778 01779 colvar%ss=r 01780 fi(:)= 2.0_dp*xi 01781 CALL put_derivative(colvar, 1, fi) 01782 01783 END SUBROUTINE xyz_diag_colvar 01784 01785 ! ***************************************************************************** 01789 SUBROUTINE xyz_outerdiag_colvar(colvar,cell,subsys,particles,error) 01790 TYPE(colvar_type), POINTER :: colvar 01791 TYPE(cell_type), POINTER :: cell 01792 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 01793 TYPE(particle_type), DIMENSION(:), 01794 OPTIONAL, POINTER :: particles 01795 TYPE(cp_error_type), INTENT(inout) :: error 01796 01797 CHARACTER(len=*), PARAMETER :: routineN = 'xyz_outerdiag_colvar', 01798 routineP = moduleN//':'//routineN 01799 01800 INTEGER :: i, k, l 01801 LOGICAL :: failure 01802 REAL(dp) :: fi(3,2), r, r0(3), ss(3), 01803 xi(3,2), xpi(3) 01804 TYPE(particle_list_type), POINTER :: particles_i 01805 TYPE(particle_type), DIMENSION(:), 01806 POINTER :: my_particles 01807 01808 failure=.FALSE. 01809 NULLIFY(particles_i) 01810 01811 CPPrecondition(colvar%type_id==xyz_outerdiag_colvar_id,cp_failure_level,routineP,error,failure) 01812 IF (PRESENT(particles)) THEN 01813 my_particles => particles 01814 ELSE 01815 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 01816 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 01817 my_particles => particles_i%els 01818 END IF 01819 DO k = 1, 2 01820 i=colvar%xyz_outerdiag_param%i_atoms(k) 01821 ! Atom coordinates 01822 CALL get_coordinates(colvar, i, xpi, my_particles) 01823 r0 = colvar%xyz_outerdiag_param%r0(:,k) 01824 IF (ALL(colvar%xyz_outerdiag_param%r0(:,k)==HUGE(0.0_dp))) r0=xpi 01825 01826 IF(colvar%xyz_outerdiag_param%use_pbc)THEN 01827 ss=MATMUL(cell%h_inv,xpi-r0) 01828 ss=ss-NINT(ss) 01829 xi(:,k)=MATMUL(cell%hmat,ss) 01830 ELSE 01831 xi(:,k)=xpi-r0 01832 END IF 01833 01834 SELECT CASE (colvar%xyz_outerdiag_param%components(k)) 01835 CASE(do_clv_x) 01836 xi(2,k)=0.0_dp 01837 xi(3,k)=0.0_dp 01838 CASE(do_clv_y) 01839 xi(1,k)=0.0_dp 01840 xi(3,k)=0.0_dp 01841 CASE(do_clv_z) 01842 xi(1,k)=0.0_dp 01843 xi(2,k)=0.0_dp 01844 CASE(do_clv_xy) 01845 xi(3,k)=0.0_dp 01846 CASE(do_clv_xz) 01847 xi(2,k)=0.0_dp 01848 CASE(do_clv_yz) 01849 xi(1,k)=0.0_dp 01850 CASE DEFAULT 01851 ! do_clv_xyz 01852 END SELECT 01853 END DO 01854 01855 r =0.0_dp 01856 fi=0.0_dp 01857 DO i =1,3 01858 DO l =1,3 01859 IF (xi(l,1)/=0.0_dp) fi(l,1)=fi(l,1)+xi(i,2) 01860 r = r + xi(l,1)*xi(i,2) 01861 END DO 01862 IF (xi(i,2)/=0.0_dp) fi(i,2)=SUM(xi(:,1)) 01863 END DO 01864 01865 colvar%ss=r 01866 CALL put_derivative(colvar, 1, fi(:,1)) 01867 CALL put_derivative(colvar, 2, fi(:,2)) 01868 01869 END SUBROUTINE xyz_outerdiag_colvar 01870 01871 ! ***************************************************************************** 01877 SUBROUTINE u_colvar(colvar,force_env,subsys,error) 01878 TYPE(colvar_type), POINTER :: colvar 01879 TYPE(force_env_type), OPTIONAL, POINTER :: force_env 01880 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 01881 TYPE(cp_error_type), INTENT(inout) :: error 01882 01883 CHARACTER(len=*), PARAMETER :: routineN = 'u_colvar', 01884 routineP = moduleN//':'//routineN 01885 01886 CHARACTER(LEN=default_path_length) :: coupling_function 01887 CHARACTER(LEN=default_string_length) :: def_error, this_error 01888 CHARACTER(LEN=default_string_length), 01889 DIMENSION(:), POINTER :: parameters 01890 INTEGER :: iatom, iforce_eval, 01891 iparticle, jparticle, natom, 01892 natom_iforce, nforce_eval, 01893 stat 01894 INTEGER, DIMENSION(:), POINTER :: glob_natoms, map_index 01895 LOGICAL :: failure 01896 REAL(dp) :: dedf, dx, err, fi(3), lerr, 01897 potential_energy 01898 REAL(KIND=dp), DIMENSION(:), POINTER :: values 01899 TYPE(cp_subsys_p_type), DIMENSION(:), 01900 POINTER :: subsystems 01901 TYPE(cp_subsys_type), POINTER :: subsys_main 01902 TYPE(mixed_force_type), DIMENSION(:), 01903 POINTER :: global_forces 01904 TYPE(particle_list_p_type), 01905 DIMENSION(:), POINTER :: particles 01906 TYPE(particle_list_type), POINTER :: particles_main 01907 TYPE(section_vals_type), POINTER :: force_env_section, 01908 mapping_section, wrk_section 01909 01910 failure = .FALSE. 01911 IF (PRESENT(force_env)) THEN 01912 NULLIFY(particles_main, subsys_main) 01913 CALL force_env_get(force_env=force_env,subsys=subsys_main, error=error) 01914 CALL cp_subsys_get(subsys=subsys_main, particles=particles_main,error=error) 01915 natom = SIZE(particles_main%els) 01916 colvar%n_atom_s = natom 01917 colvar%u_param%natom = natom 01918 CALL reallocate(colvar%i_atom, 1,natom) 01919 CALL reallocate(colvar%dsdr,1,3,1,natom) 01920 DO iatom=1,natom 01921 colvar%i_atom(iatom)=iatom 01922 ENDDO 01923 01924 IF (.NOT.ASSOCIATED(colvar%u_param%mixed_energy_section)) THEN 01925 CALL force_env_get(force_env,potential_energy=potential_energy,error=error) 01926 colvar%ss = potential_energy 01927 01928 DO iatom=1,natom 01929 ! store derivative 01930 fi(:)=-particles_main%els(iatom)%f 01931 CALL put_derivative(colvar, iatom, fi) 01932 ENDDO 01933 ELSE 01934 CALL cp_assert(force_env%in_use==use_mixed_force,cp_fatal_level,-300,routineP,& 01935 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)//& 01936 ' A combination of mixed force_eval energies has been requested as '//& 01937 ' collective variable, but the MIXED env is not in use! Aborting.',& 01938 error=error,only_ionode=.TRUE.) 01939 CALL force_env_get(force_env, force_env_section=force_env_section, error=error) 01940 mapping_section => section_vals_get_subs_vals(force_env_section,"MIXED%MAPPING",error=error) 01941 NULLIFY(values, parameters, subsystems, particles, global_forces, map_index, glob_natoms) 01942 nforce_eval = SIZE(force_env%sub_force_env) 01943 ALLOCATE(glob_natoms(nforce_eval), stat=stat) 01944 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01945 ALLOCATE(subsystems(nforce_eval), stat=stat) 01946 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01947 ALLOCATE(particles(nforce_eval), stat=stat) 01948 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01949 ! Local Info to sync 01950 ALLOCATE(global_forces(nforce_eval), stat=stat) 01951 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01952 01953 glob_natoms=0 01954 DO iforce_eval = 1, nforce_eval 01955 NULLIFY(subsystems(iforce_eval)%subsys, particles(iforce_eval)%list) 01956 IF (.NOT.ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE 01957 ! Get all available subsys 01958 CALL force_env_get(force_env=force_env%sub_force_env(iforce_eval)%force_env,& 01959 subsys=subsystems(iforce_eval)%subsys,error=error) 01960 ! Get available particles 01961 CALL cp_subsys_get(subsys=subsystems(iforce_eval)%subsys,& 01962 particles=particles(iforce_eval)%list,error=error) 01963 01964 ! Get Mapping index array 01965 natom_iforce = SIZE(particles(iforce_eval)%list%els) 01966 01967 ! Only the rank 0 process collect info for each computation 01968 IF ( force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos==& 01969 force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN 01970 glob_natoms(iforce_eval) = natom_iforce 01971 END IF 01972 END DO 01973 01974 ! Handling Parallel execution 01975 CALL mp_sync(force_env%para_env%group) 01976 CALL mp_sum(glob_natoms, force_env%para_env%group) 01977 01978 ! Transfer forces 01979 DO iforce_eval = 1, nforce_eval 01980 ALLOCATE(global_forces(iforce_eval)%forces(3,glob_natoms(iforce_eval)),stat=stat) 01981 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01982 global_forces(iforce_eval)%forces = 0.0_dp 01983 IF (ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) THEN 01984 IF ( force_env%sub_force_env(iforce_eval)%force_env%para_env%mepos==& 01985 force_env%sub_force_env(iforce_eval)%force_env%para_env%source) THEN 01986 ! Forces 01987 DO iparticle = 1, glob_natoms(iforce_eval) 01988 global_forces(iforce_eval)%forces(:,iparticle) = & 01989 particles(iforce_eval)%list%els(iparticle)%f 01990 END DO 01991 END IF 01992 END IF 01993 CALL mp_sum(global_forces(iforce_eval)%forces, force_env%para_env%group) 01994 END DO 01995 01996 wrk_section => colvar%u_param%mixed_energy_section 01997 ! Support any number of force_eval sections 01998 CALL get_generic_info(wrk_section, "ENERGY_FUNCTION", coupling_function, parameters,& 01999 values, force_env%mixed_env%energies, error=error) 02000 CALL initf(1) 02001 CALL parsef(1,TRIM(coupling_function),parameters) 02002 ! Store the value of the COLVAR 02003 colvar%ss = evalf(1,values) 02004 CPPrecondition(EvalErrType<=0,cp_failure_level,routineP,error,failure) 02005 02006 DO iforce_eval = 1, nforce_eval 02007 CALL section_vals_val_get(wrk_section,"DX",r_val=dx,error=error) 02008 CALL section_vals_val_get(wrk_section,"ERROR_LIMIT",r_val=lerr,error=error) 02009 dedf = evalfd(1,iforce_eval,values,dx,err) 02010 IF (ABS(err)>lerr) THEN 02011 WRITE(this_error,"(A,G12.6,A)")"(",err,")" 02012 WRITE(def_error,"(A,G12.6,A)")"(",lerr,")" 02013 CALL compress(this_error,.TRUE.) 02014 CALL compress(def_error,.TRUE.) 02015 CALL cp_assert(.FALSE.,cp_warning_level,-300,routineP,& 02016 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)//& 02017 ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'//& 02018 TRIM(def_error)//' .',error=error,only_ionode=.TRUE.) 02019 END IF 02020 ! General Mapping of forces... 02021 ! First: Get Mapping index array 02022 CALL get_subsys_map_index(mapping_section, glob_natoms(iforce_eval), iforce_eval,& 02023 nforce_eval, map_index, error) 02024 02025 ! Second: store derivatives 02026 DO iparticle = 1, glob_natoms(iforce_eval) 02027 jparticle = map_index(iparticle) 02028 fi= - dedf * global_forces(iforce_eval)%forces(:,iparticle) 02029 CALL put_derivative(colvar, jparticle, fi) 02030 END DO 02031 ! Deallocate map_index array 02032 IF (ASSOCIATED(map_index)) THEN 02033 DEALLOCATE(map_index, stat=stat) 02034 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02035 END IF 02036 END DO 02037 CALL finalizef() 02038 DO iforce_eval = 1, nforce_eval 02039 DEALLOCATE(global_forces(iforce_eval)%forces,stat=stat) 02040 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02041 END DO 02042 DEALLOCATE(glob_natoms, stat=stat) 02043 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02044 DEALLOCATE(values, stat=stat) 02045 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02046 DEALLOCATE(parameters, stat=stat) 02047 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02048 DEALLOCATE(global_forces, stat=stat) 02049 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02050 DEALLOCATE(subsystems, stat=stat) 02051 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02052 DEALLOCATE(particles, stat=stat) 02053 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 02054 END IF 02055 ELSE 02056 CALL cp_unimplemented_error(fromWhere=routineP, & 02057 message="need force_env!", error=error, error_level=cp_failure_level) 02058 ENDIF 02059 END SUBROUTINE u_colvar 02060 02061 ! ***************************************************************************** 02065 SUBROUTINE plane_distance_colvar(colvar,cell,subsys,particles,error) 02066 02067 TYPE(colvar_type), POINTER :: colvar 02068 TYPE(cell_type), POINTER :: cell 02069 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02070 TYPE(particle_type), DIMENSION(:), 02071 OPTIONAL, POINTER :: particles 02072 TYPE(cp_error_type), INTENT(inout) :: error 02073 02074 CHARACTER(len=*), PARAMETER :: routineN = 'plane_distance_colvar', 02075 routineP = moduleN//':'//routineN 02076 02077 INTEGER :: i, j, k, l 02078 LOGICAL :: failure 02079 REAL(dp) :: a, b, dsdxpn(3), dxpndxi(3,3), dxpndxj(3,3), dxpndxk(3,3), 02080 fi(3), fj(3), fk(3), fl(3), r12, ri(3), rj(3), rk(3), rl(3), ss(3), 02081 xpij(3), xpkj(3), xpl(3), xpn(3) 02082 TYPE(particle_list_type), POINTER :: particles_i 02083 TYPE(particle_type), DIMENSION(:), 02084 POINTER :: my_particles 02085 02086 failure=.FALSE. 02087 NULLIFY(particles_i) 02088 02089 CPPrecondition(colvar%type_id==plane_distance_colvar_id,cp_failure_level,routineP,error,failure) 02090 IF (PRESENT(particles)) THEN 02091 my_particles => particles 02092 ELSE 02093 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02094 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02095 my_particles => particles_i%els 02096 END IF 02097 i=colvar%plane_distance_param%plane(1) 02098 j=colvar%plane_distance_param%plane(2) 02099 k=colvar%plane_distance_param%plane(3) 02100 l=colvar%plane_distance_param%point 02101 ! Get coordinates of atoms or points 02102 CALL get_coordinates(colvar, i, ri, my_particles) 02103 CALL get_coordinates(colvar, j, rj, my_particles) 02104 CALL get_coordinates(colvar, k, rk, my_particles) 02105 CALL get_coordinates(colvar, l, rl, my_particles) 02106 xpij=ri-rj 02107 xpkj=rk-rj 02108 xpl=rl-(ri+rj+rk)/3.0_dp 02109 IF (colvar%plane_distance_param%use_pbc) THEN 02110 ! xpij 02111 ss=MATMUL(cell%h_inv,ri-rj) 02112 ss=ss-NINT(ss) 02113 xpij=MATMUL(cell%hmat,ss) 02114 ! xpkj 02115 ss=MATMUL(cell%h_inv,rk-rj) 02116 ss=ss-NINT(ss) 02117 xpkj=MATMUL(cell%hmat,ss) 02118 ! xpl 02119 ss=MATMUL(cell%h_inv,rl-(ri+rj+rk)/3.0_dp) 02120 ss=ss-NINT(ss) 02121 xpl=MATMUL(cell%hmat,ss) 02122 END IF 02123 ! xpn 02124 xpn(1) = xpij(2)*xpkj(3)-xpij(3)*xpkj(2) 02125 xpn(2) = xpij(3)*xpkj(1)-xpij(1)*xpkj(3) 02126 xpn(3) = xpij(1)*xpkj(2)-xpij(2)*xpkj(1) 02127 a = DOT_PRODUCT(xpn,xpn) 02128 b = DOT_PRODUCT(xpl,xpn) 02129 r12=SQRT(a) 02130 colvar%ss=b/r12 02131 dsdxpn(1) = xpl(1)/r12 - b*xpn(1)/(r12*a) 02132 dsdxpn(2) = xpl(2)/r12 - b*xpn(2)/(r12*a) 02133 dsdxpn(3) = xpl(3)/r12 - b*xpn(3)/(r12*a) 02134 ! 02135 dxpndxi(1,1)= 0.0_dp 02136 dxpndxi(1,2)= 1.0_dp * xpkj(3) 02137 dxpndxi(1,3)= -1.0_dp * xpkj(2) 02138 dxpndxi(2,1)= -1.0_dp * xpkj(3) 02139 dxpndxi(2,2)= 0.0_dp 02140 dxpndxi(2,3)= 1.0_dp * xpkj(1) 02141 dxpndxi(3,1)= 1.0_dp * xpkj(2) 02142 dxpndxi(3,2)= -1.0_dp * xpkj(1) 02143 dxpndxi(3,3)= 0.0_dp 02144 ! 02145 dxpndxj(1,1)= 0.0_dp 02146 dxpndxj(1,2)= -1.0_dp * xpkj(3) + xpij(3) 02147 dxpndxj(1,3)= -1.0_dp * xpij(2) + xpkj(2) 02148 dxpndxj(2,1)= -1.0_dp * xpij(3) + xpkj(3) 02149 dxpndxj(2,2)= 0.0_dp 02150 dxpndxj(2,3)= -1.0_dp * xpkj(1) + xpij(1) 02151 dxpndxj(3,1)= -1.0_dp * xpkj(2) + xpij(2) 02152 dxpndxj(3,2)= -1.0_dp * xpij(1) + xpkj(1) 02153 dxpndxj(3,3)= 0.0_dp 02154 ! 02155 dxpndxk(1,1)= 0.0_dp 02156 dxpndxk(1,2)= -1.0_dp * xpij(3) 02157 dxpndxk(1,3)= 1.0_dp * xpij(2) 02158 dxpndxk(2,1)= 1.0_dp * xpij(3) 02159 dxpndxk(2,2)= 0.0_dp 02160 dxpndxk(2,3)= -1.0_dp * xpij(1) 02161 dxpndxk(3,1)= -1.0_dp * xpij(2) 02162 dxpndxk(3,2)= 1.0_dp * xpij(1) 02163 dxpndxk(3,3)= 0.0_dp 02164 ! 02165 fi(:)=MATMUL(dsdxpn,dxpndxi)-xpn/(3.0_dp*r12) 02166 fj(:)=MATMUL(dsdxpn,dxpndxj)-xpn/(3.0_dp*r12) 02167 fk(:)=MATMUL(dsdxpn,dxpndxk)-xpn/(3.0_dp*r12) 02168 fl(:)=xpn/r12 02169 ! Transfer derivatives on atoms 02170 CALL put_derivative(colvar, 1, fi) 02171 CALL put_derivative(colvar, 2, fj) 02172 CALL put_derivative(colvar, 3, fk) 02173 CALL put_derivative(colvar, 4, fl) 02174 02175 END SUBROUTINE plane_distance_colvar 02176 02177 02178 ! ***************************************************************************** 02183 SUBROUTINE plane_plane_angle_colvar(colvar,cell,subsys,particles,error) 02184 02185 TYPE(colvar_type), POINTER :: colvar 02186 TYPE(cell_type), POINTER :: cell 02187 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02188 TYPE(particle_type), DIMENSION(:), 02189 OPTIONAL, POINTER :: particles 02190 TYPE(cp_error_type), INTENT(inout) :: error 02191 02192 CHARACTER(len=*), PARAMETER :: routineN = 'plane_plane_angle_colvar', 02193 routineP = moduleN//':'//routineN 02194 02195 INTEGER :: i1, i2, j1, j2, k1, k2, np 02196 LOGICAL :: check, failure 02197 REAL(dp) :: a1, a2, d, dnorm_dxpn(3), dprod12_dxpn(3), dsdxpn(3), 02198 dt_dxpn(3), dxpndxi(3,3), dxpndxj(3,3), dxpndxk(3,3), fi(3), fj(3), 02199 fk(3), fmod, norm1, norm2, prod_12, ri1(3), ri2(3), rj1(3), rj2(3), 02200 rk1(3), rk2(3), ss(3), t, xpij1(3), xpij2(3), xpkj1(3), xpkj2(3), 02201 xpn1(3), xpn2(3) 02202 TYPE(particle_list_type), POINTER :: particles_i 02203 TYPE(particle_type), DIMENSION(:), 02204 POINTER :: my_particles 02205 02206 failure=.FALSE. 02207 NULLIFY(particles_i) 02208 02209 check = colvar%type_id==plane_plane_angle_colvar_id 02210 CPPrecondition(check,cp_failure_level,routineP,error,failure) 02211 IF (PRESENT(particles)) THEN 02212 my_particles => particles 02213 ELSE 02214 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02215 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02216 my_particles => particles_i%els 02217 END IF 02218 02219 ! Plane 1 02220 IF (colvar%plane_plane_angle_param%plane1%type_of_def==plane_def_atoms) THEN 02221 i1=colvar%plane_plane_angle_param%plane1%points(1) 02222 j1=colvar%plane_plane_angle_param%plane1%points(2) 02223 k1=colvar%plane_plane_angle_param%plane1%points(3) 02224 02225 ! Get coordinates of atoms or points 02226 CALL get_coordinates(colvar, i1, ri1, my_particles) 02227 CALL get_coordinates(colvar, j1, rj1, my_particles) 02228 CALL get_coordinates(colvar, k1, rk1, my_particles) 02229 02230 ! xpij 02231 ss=MATMUL(cell%h_inv,ri1-rj1) 02232 ss=ss-NINT(ss) 02233 xpij1=MATMUL(cell%hmat,ss) 02234 02235 ! xpkj 02236 ss=MATMUL(cell%h_inv,rk1-rj1) 02237 ss=ss-NINT(ss) 02238 xpkj1=MATMUL(cell%hmat,ss) 02239 02240 ! xpn 02241 xpn1(1) = xpij1(2)*xpkj1(3)-xpij1(3)*xpkj1(2) 02242 xpn1(2) = xpij1(3)*xpkj1(1)-xpij1(1)*xpkj1(3) 02243 xpn1(3) = xpij1(1)*xpkj1(2)-xpij1(2)*xpkj1(1) 02244 ELSE 02245 xpn1 = colvar%plane_plane_angle_param%plane1%normal_vec 02246 END IF 02247 a1 = DOT_PRODUCT(xpn1,xpn1) 02248 norm1 = SQRT(a1) 02249 CPPrecondition(norm1/=0.0_dp,cp_failure_level,routineP,error,failure) 02250 02251 ! Plane 2 02252 IF (colvar%plane_plane_angle_param%plane2%type_of_def==plane_def_atoms) THEN 02253 i2=colvar%plane_plane_angle_param%plane2%points(1) 02254 j2=colvar%plane_plane_angle_param%plane2%points(2) 02255 k2=colvar%plane_plane_angle_param%plane2%points(3) 02256 02257 ! Get coordinates of atoms or points 02258 CALL get_coordinates(colvar, i2, ri2, my_particles) 02259 CALL get_coordinates(colvar, j2, rj2, my_particles) 02260 CALL get_coordinates(colvar, k2, rk2, my_particles) 02261 02262 ! xpij 02263 ss=MATMUL(cell%h_inv,ri2-rj2) 02264 ss=ss-NINT(ss) 02265 xpij2=MATMUL(cell%hmat,ss) 02266 02267 ! xpkj 02268 ss=MATMUL(cell%h_inv,rk2-rj2) 02269 ss=ss-NINT(ss) 02270 xpkj2=MATMUL(cell%hmat,ss) 02271 02272 ! xpn 02273 xpn2(1) = xpij2(2)*xpkj2(3)-xpij2(3)*xpkj2(2) 02274 xpn2(2) = xpij2(3)*xpkj2(1)-xpij2(1)*xpkj2(3) 02275 xpn2(3) = xpij2(1)*xpkj2(2)-xpij2(2)*xpkj2(1) 02276 ELSE 02277 xpn2 = colvar%plane_plane_angle_param%plane2%normal_vec 02278 END IF 02279 a2 = DOT_PRODUCT(xpn2,xpn2) 02280 norm2 = SQRT(a2) 02281 CPPrecondition(norm2/=0.0_dp,cp_failure_level,routineP,error,failure) 02282 02283 ! The value of the angle is defined only between 0 and Pi 02284 prod_12 = DOT_PRODUCT(xpn1,xpn2) 02285 02286 d = norm1*norm2 02287 t = prod_12/d 02288 t = MIN(1.0_dp,ABS(t))*SIGN(1.0_dp,t) 02289 colvar%ss= ACOS(t) 02290 02291 IF ((ABS(colvar%ss).LT.tolerance_acos).OR.(ABS(colvar%ss-pi).LT.tolerance_acos)) THEN 02292 fmod = 0.0_dp 02293 ELSE 02294 fmod = - 1.0_dp / SIN(colvar%ss) 02295 ENDIF 02296 ! Compute derivatives 02297 np = 0 02298 ! Plane 1 02299 IF (colvar%plane_plane_angle_param%plane1%type_of_def==plane_def_atoms) THEN 02300 dprod12_dxpn = xpn2 02301 dnorm_dxpn = 1.0_dp/norm1*xpn1 02302 dt_dxpn = (dprod12_dxpn * d - prod_12 * dnorm_dxpn * norm2)/d**2 02303 02304 dsdxpn(1) = fmod * dt_dxpn(1) 02305 dsdxpn(2) = fmod * dt_dxpn(2) 02306 dsdxpn(3) = fmod * dt_dxpn(3) 02307 ! 02308 dxpndxi(1,1)= 0.0_dp 02309 dxpndxi(1,2)= 1.0_dp * xpkj1(3) 02310 dxpndxi(1,3)= -1.0_dp * xpkj1(2) 02311 dxpndxi(2,1)= -1.0_dp * xpkj1(3) 02312 dxpndxi(2,2)= 0.0_dp 02313 dxpndxi(2,3)= 1.0_dp * xpkj1(1) 02314 dxpndxi(3,1)= 1.0_dp * xpkj1(2) 02315 dxpndxi(3,2)= -1.0_dp * xpkj1(1) 02316 dxpndxi(3,3)= 0.0_dp 02317 ! 02318 dxpndxj(1,1)= 0.0_dp 02319 dxpndxj(1,2)= -1.0_dp * xpkj1(3) + xpij1(3) 02320 dxpndxj(1,3)= -1.0_dp * xpij1(2) + xpkj1(2) 02321 dxpndxj(2,1)= -1.0_dp * xpij1(3) + xpkj1(3) 02322 dxpndxj(2,2)= 0.0_dp 02323 dxpndxj(2,3)= -1.0_dp * xpkj1(1) + xpij1(1) 02324 dxpndxj(3,1)= -1.0_dp * xpkj1(2) + xpij1(2) 02325 dxpndxj(3,2)= -1.0_dp * xpij1(1) + xpkj1(1) 02326 dxpndxj(3,3)= 0.0_dp 02327 ! 02328 dxpndxk(1,1)= 0.0_dp 02329 dxpndxk(1,2)= -1.0_dp * xpij1(3) 02330 dxpndxk(1,3)= 1.0_dp * xpij1(2) 02331 dxpndxk(2,1)= 1.0_dp * xpij1(3) 02332 dxpndxk(2,2)= 0.0_dp 02333 dxpndxk(2,3)= -1.0_dp * xpij1(1) 02334 dxpndxk(3,1)= -1.0_dp * xpij1(2) 02335 dxpndxk(3,2)= 1.0_dp * xpij1(1) 02336 dxpndxk(3,3)= 0.0_dp 02337 ! 02338 fi=MATMUL(dsdxpn,dxpndxi) 02339 fj=MATMUL(dsdxpn,dxpndxj) 02340 fk=MATMUL(dsdxpn,dxpndxk) 02341 02342 ! Transfer derivatives on atoms 02343 CALL put_derivative(colvar, np+1, fi) 02344 CALL put_derivative(colvar, np+2, fj) 02345 CALL put_derivative(colvar, np+3, fk) 02346 np = 3 02347 END IF 02348 02349 ! Plane 2 02350 IF (colvar%plane_plane_angle_param%plane2%type_of_def==plane_def_atoms) THEN 02351 dprod12_dxpn = xpn1 02352 dnorm_dxpn = 1.0_dp/norm2*xpn2 02353 dt_dxpn = (dprod12_dxpn * d - prod_12 * dnorm_dxpn * norm1)/d**2 02354 02355 dsdxpn(1) = fmod * dt_dxpn(1) 02356 dsdxpn(2) = fmod * dt_dxpn(2) 02357 dsdxpn(3) = fmod * dt_dxpn(3) 02358 ! 02359 dxpndxi(1,1)= 0.0_dp 02360 dxpndxi(1,2)= 1.0_dp * xpkj1(3) 02361 dxpndxi(1,3)= -1.0_dp * xpkj1(2) 02362 dxpndxi(2,1)= -1.0_dp * xpkj1(3) 02363 dxpndxi(2,2)= 0.0_dp 02364 dxpndxi(2,3)= 1.0_dp * xpkj1(1) 02365 dxpndxi(3,1)= 1.0_dp * xpkj1(2) 02366 dxpndxi(3,2)= -1.0_dp * xpkj1(1) 02367 dxpndxi(3,3)= 0.0_dp 02368 ! 02369 dxpndxj(1,1)= 0.0_dp 02370 dxpndxj(1,2)= -1.0_dp * xpkj1(3) + xpij1(3) 02371 dxpndxj(1,3)= -1.0_dp * xpij1(2) + xpkj1(2) 02372 dxpndxj(2,1)= -1.0_dp * xpij1(3) + xpkj1(3) 02373 dxpndxj(2,2)= 0.0_dp 02374 dxpndxj(2,3)= -1.0_dp * xpkj1(1) + xpij1(1) 02375 dxpndxj(3,1)= -1.0_dp * xpkj1(2) + xpij1(2) 02376 dxpndxj(3,2)= -1.0_dp * xpij1(1) + xpkj1(1) 02377 dxpndxj(3,3)= 0.0_dp 02378 ! 02379 dxpndxk(1,1)= 0.0_dp 02380 dxpndxk(1,2)= -1.0_dp * xpij1(3) 02381 dxpndxk(1,3)= 1.0_dp * xpij1(2) 02382 dxpndxk(2,1)= 1.0_dp * xpij1(3) 02383 dxpndxk(2,2)= 0.0_dp 02384 dxpndxk(2,3)= -1.0_dp * xpij1(1) 02385 dxpndxk(3,1)= -1.0_dp * xpij1(2) 02386 dxpndxk(3,2)= 1.0_dp * xpij1(1) 02387 dxpndxk(3,3)= 0.0_dp 02388 ! 02389 fi=MATMUL(dsdxpn,dxpndxi) 02390 fj=MATMUL(dsdxpn,dxpndxj) 02391 fk=MATMUL(dsdxpn,dxpndxk) 02392 02393 ! Transfer derivatives on atoms 02394 CALL put_derivative(colvar, np+1, fi) 02395 CALL put_derivative(colvar, np+2, fj) 02396 CALL put_derivative(colvar, np+3, fk) 02397 END IF 02398 02399 END SUBROUTINE plane_plane_angle_colvar 02400 02401 ! ***************************************************************************** 02405 SUBROUTINE rotation_colvar(colvar,cell,subsys,particles,error) 02406 TYPE(colvar_type), POINTER :: colvar 02407 TYPE(cell_type), POINTER :: cell 02408 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02409 TYPE(particle_type), DIMENSION(:), 02410 OPTIONAL, POINTER :: particles 02411 TYPE(cp_error_type), INTENT(inout) :: error 02412 02413 CHARACTER(len=*), PARAMETER :: routineN = 'rotation_colvar', 02414 routineP = moduleN//':'//routineN 02415 02416 INTEGER :: i, idum 02417 LOGICAL :: failure 02418 REAL(dp) :: a, b, fmod, t0, t1, t2, t3, 02419 xdum(3), xij(3), xkj(3) 02420 REAL(KIND=dp) :: dp1b1(3), dp1b2(3), dp2b1(3), 02421 dp2b2(3), ss(3), xp1b1(3), 02422 xp1b2(3), xp2b1(3), xp2b2(3) 02423 TYPE(particle_list_type), POINTER :: particles_i 02424 TYPE(particle_type), DIMENSION(:), 02425 POINTER :: my_particles 02426 02427 failure=.FALSE. 02428 NULLIFY(particles_i) 02429 02430 CPPrecondition(colvar%type_id==rotation_colvar_id,cp_failure_level,routineP,error,failure) 02431 IF (PRESENT(particles)) THEN 02432 my_particles => particles 02433 ELSE 02434 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02435 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02436 my_particles => particles_i%els 02437 END IF 02438 i = colvar%rotation_param%i_at1_bond1 02439 CALL get_coordinates(colvar, i, xp1b1, my_particles) 02440 i = colvar%rotation_param%i_at2_bond1 02441 CALL get_coordinates(colvar, i, xp2b1, my_particles) 02442 i = colvar%rotation_param%i_at1_bond2 02443 CALL get_coordinates(colvar, i, xp1b2, my_particles) 02444 i = colvar%rotation_param%i_at2_bond2 02445 CALL get_coordinates(colvar, i, xp2b2, my_particles) 02446 ! xij 02447 ss=MATMUL(cell%h_inv,xp1b1-xp2b1) 02448 ss=ss-NINT(ss) 02449 xij=MATMUL(cell%hmat,ss) 02450 ! xkj 02451 ss=MATMUL(cell%h_inv,xp1b2-xp2b2) 02452 ss=ss-NINT(ss) 02453 xkj=MATMUL(cell%hmat,ss) 02454 ! evaluation of the angle.. 02455 a = SQRT(DOT_PRODUCT(xij,xij)) 02456 b = SQRT(DOT_PRODUCT(xkj,xkj)) 02457 t0 = 1.0_dp/(a*b) 02458 t1 = 1.0_dp/(a**3.0_dp*b) 02459 t2 = 1.0_dp/(a*b**3.0_dp) 02460 t3 = DOT_PRODUCT(xij,xkj) 02461 colvar%ss = ACOS(t3*t0) 02462 IF ((ABS(colvar%ss).LT.tolerance_acos).OR.(ABS(colvar%ss-pi).LT.tolerance_acos)) THEN 02463 fmod = 0.0_dp 02464 ELSE 02465 fmod = - 1.0_dp / SIN(colvar%ss) 02466 ENDIF 02467 dp1b1 = xkj(:)*t0-xij(:)*t1*t3 02468 dp2b1 = -xkj(:)*t0+xij(:)*t1*t3 02469 dp1b2 = xij(:)*t0-xkj(:)*t2*t3 02470 dp2b2 = -xij(:)*t0+xkj(:)*t2*t3 02471 02472 xdum = dp1b1 * fmod 02473 idum = colvar%rotation_param%i_at1_bond1 02474 CALL put_derivative(colvar, idum, xdum) 02475 xdum = dp2b1 * fmod 02476 idum = colvar%rotation_param%i_at2_bond1 02477 CALL put_derivative(colvar, idum, xdum) 02478 xdum = dp1b2 * fmod 02479 idum = colvar%rotation_param%i_at1_bond2 02480 CALL put_derivative(colvar, idum, xdum) 02481 xdum = dp2b2 * fmod 02482 idum = colvar%rotation_param%i_at2_bond2 02483 CALL put_derivative(colvar, idum, xdum) 02484 02485 END SUBROUTINE rotation_colvar 02486 02487 ! ***************************************************************************** 02492 SUBROUTINE dfunct_colvar(colvar,cell,subsys,particles,error) 02493 TYPE(colvar_type), POINTER :: colvar 02494 TYPE(cell_type), POINTER :: cell 02495 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02496 TYPE(particle_type), DIMENSION(:), 02497 OPTIONAL, POINTER :: particles 02498 TYPE(cp_error_type), INTENT(inout) :: error 02499 02500 CHARACTER(len=*), PARAMETER :: routineN = 'dfunct_colvar', 02501 routineP = moduleN//':'//routineN 02502 02503 INTEGER :: i, j, k, l 02504 LOGICAL :: failure 02505 REAL(dp) :: fi(3), fj(3), fk(3), fl(3), 02506 r12, r34, ss(3), xij(3), 02507 xkl(3), xpi(3), xpj(3), 02508 xpk(3), xpl(3) 02509 TYPE(particle_list_type), POINTER :: particles_i 02510 TYPE(particle_type), DIMENSION(:), 02511 POINTER :: my_particles 02512 02513 failure=.FALSE. 02514 NULLIFY(particles_i) 02515 02516 CPPrecondition(colvar%type_id==dfunct_colvar_id,cp_failure_level,routineP,error,failure) 02517 IF (PRESENT(particles)) THEN 02518 my_particles => particles 02519 ELSE 02520 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02521 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02522 my_particles => particles_i%els 02523 END IF 02524 i=colvar%dfunct_param%i_at_dfunct(1) 02525 j=colvar%dfunct_param%i_at_dfunct(2) 02526 ! First bond 02527 CALL get_coordinates(colvar, i, xpi, my_particles) 02528 CALL get_coordinates(colvar, j, xpj, my_particles) 02529 IF(colvar%dfunct_param%use_pbc)THEN 02530 ss=MATMUL(cell%h_inv,xpi-xpj) 02531 ss=ss-NINT(ss) 02532 xij=MATMUL(cell%hmat,ss) 02533 ELSE 02534 xij=xpi-xpj 02535 END IF 02536 r12=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 02537 ! Second bond 02538 k=colvar%dfunct_param%i_at_dfunct(3) 02539 l=colvar%dfunct_param%i_at_dfunct(4) 02540 CALL get_coordinates(colvar, k, xpk, my_particles) 02541 CALL get_coordinates(colvar, l, xpl, my_particles) 02542 IF(colvar%dfunct_param%use_pbc)THEN 02543 ss=MATMUL(cell%h_inv,xpk-xpl) 02544 ss=ss-NINT(ss) 02545 xkl=MATMUL(cell%hmat,ss) 02546 ELSE 02547 xkl=xpk-xpl 02548 END IF 02549 r34=SQRT(xkl(1)**2+xkl(2)**2+xkl(3)**2) 02550 ! 02551 colvar%ss=r12+colvar%dfunct_param%coeff*r34 02552 fi(:)= xij/r12 02553 fj(:)=-xij/r12 02554 fk(:)=colvar%dfunct_param%coeff*xkl/r34 02555 fl(:)=-colvar%dfunct_param%coeff*xkl/r34 02556 CALL put_derivative(colvar, 1, fi) 02557 CALL put_derivative(colvar, 2, fj) 02558 CALL put_derivative(colvar, 3, fk) 02559 CALL put_derivative(colvar, 4, fl) 02560 02561 END SUBROUTINE dfunct_colvar 02562 02563 ! ***************************************************************************** 02567 SUBROUTINE angle_colvar(colvar,cell,subsys,particles,error) 02568 TYPE(colvar_type), POINTER :: colvar 02569 TYPE(cell_type), POINTER :: cell 02570 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02571 TYPE(particle_type), DIMENSION(:), 02572 OPTIONAL, POINTER :: particles 02573 TYPE(cp_error_type), INTENT(inout) :: error 02574 02575 CHARACTER(len=*), PARAMETER :: routineN = 'angle_colvar', 02576 routineP = moduleN//':'//routineN 02577 02578 INTEGER :: i, j, k 02579 LOGICAL :: failure 02580 REAL(dp) :: a, b, fi(3), fj(3), fk(3), 02581 fmod, ri(3), rj(3), rk(3), 02582 ss(3), t0, t1, t2, t3, 02583 xij(3), xkj(3) 02584 TYPE(particle_list_type), POINTER :: particles_i 02585 TYPE(particle_type), DIMENSION(:), 02586 POINTER :: my_particles 02587 02588 failure=.FALSE. 02589 NULLIFY(particles_i) 02590 02591 CPPrecondition(colvar%type_id==angle_colvar_id,cp_failure_level,routineP,error,failure) 02592 IF (PRESENT(particles)) THEN 02593 my_particles => particles 02594 ELSE 02595 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02596 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02597 my_particles => particles_i%els 02598 END IF 02599 i=colvar%angle_param%i_at_angle(1) 02600 j=colvar%angle_param%i_at_angle(2) 02601 k=colvar%angle_param%i_at_angle(3) 02602 CALL get_coordinates(colvar, i, ri, my_particles) 02603 CALL get_coordinates(colvar, j, rj, my_particles) 02604 CALL get_coordinates(colvar, k, rk, my_particles) 02605 ! xij 02606 ss=MATMUL(cell%h_inv,ri-rj) 02607 ss=ss-NINT(ss) 02608 xij=MATMUL(cell%hmat,ss) 02609 ! xkj 02610 ss=MATMUL(cell%h_inv,rk-rj) 02611 ss=ss-NINT(ss) 02612 xkj=MATMUL(cell%hmat,ss) 02613 ! Evaluation of the angle.. 02614 a = SQRT(DOT_PRODUCT(xij,xij)) 02615 b = SQRT(DOT_PRODUCT(xkj,xkj)) 02616 t0 = 1.0_dp/(a*b) 02617 t1 = 1.0_dp/(a**3.0_dp*b) 02618 t2 = 1.0_dp/(a*b**3.0_dp) 02619 t3 = DOT_PRODUCT(xij,xkj) 02620 colvar%ss = ACOS(t3*t0) 02621 IF ((ABS(colvar%ss).LT.tolerance_acos).OR.(ABS(colvar%ss-pi).LT.tolerance_acos)) THEN 02622 fmod = 0.0_dp 02623 ELSE 02624 fmod = - 1.0_dp / SIN(colvar%ss) 02625 ENDIF 02626 fi(:) = xkj(:)*t0-xij(:)*t1*t3 02627 fj(:) =-xkj(:)*t0+xij(:)*t1*t3-xij(:)*t0+xkj(:)*t2*t3 02628 fk(:) = xij(:)*t0-xkj(:)*t2*t3 02629 fi = fi * fmod 02630 fj = fj * fmod 02631 fk = fk * fmod 02632 CALL put_derivative(colvar, 1, fi) 02633 CALL put_derivative(colvar, 2, fj) 02634 CALL put_derivative(colvar, 3, fk) 02635 02636 END SUBROUTINE angle_colvar 02637 02638 ! ***************************************************************************** 02642 SUBROUTINE dist_colvar(colvar,cell,subsys,particles,error) 02643 TYPE(colvar_type), POINTER :: colvar 02644 TYPE(cell_type), POINTER :: cell 02645 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02646 TYPE(particle_type), DIMENSION(:), 02647 OPTIONAL, POINTER :: particles 02648 TYPE(cp_error_type), INTENT(inout) :: error 02649 02650 CHARACTER(len=*), PARAMETER :: routineN = 'dist_colvar', 02651 routineP = moduleN//':'//routineN 02652 02653 INTEGER :: i, j 02654 LOGICAL :: failure 02655 REAL(dp) :: fi(3), fj(3), r12, ss(3), 02656 xij(3), xpi(3), xpj(3) 02657 TYPE(particle_list_type), POINTER :: particles_i 02658 TYPE(particle_type), DIMENSION(:), 02659 POINTER :: my_particles 02660 02661 failure=.FALSE. 02662 NULLIFY(particles_i) 02663 02664 CPPrecondition(colvar%type_id==dist_colvar_id,cp_failure_level,routineP,error,failure) 02665 IF (PRESENT(particles)) THEN 02666 my_particles => particles 02667 ELSE 02668 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02669 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02670 my_particles => particles_i%els 02671 END IF 02672 i=colvar%dist_param%i_at 02673 j=colvar%dist_param%j_at 02674 CALL get_coordinates(colvar, i, xpi, my_particles) 02675 CALL get_coordinates(colvar, j, xpj, my_particles) 02676 ss=MATMUL(cell%h_inv,xpi-xpj) 02677 ss=ss-NINT(ss) 02678 xij=MATMUL(cell%hmat,ss) 02679 SELECT CASE (colvar%dist_param%axis_id) 02680 CASE(do_clv_x) 02681 xij(2)=0.0_dp 02682 xij(3)=0.0_dp 02683 CASE(do_clv_y) 02684 xij(1)=0.0_dp 02685 xij(3)=0.0_dp 02686 CASE(do_clv_z) 02687 xij(1)=0.0_dp 02688 xij(2)=0.0_dp 02689 CASE(do_clv_xy) 02690 xij(3)=0.0_dp 02691 CASE(do_clv_xz) 02692 xij(2)=0.0_dp 02693 CASE(do_clv_yz) 02694 xij(1)=0.0_dp 02695 CASE DEFAULT 02696 !do_clv_xyz 02697 END SELECT 02698 r12=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 02699 02700 colvar%ss=r12 02701 fi(:)=xij/r12 02702 fj(:)=-xij/r12 02703 02704 CALL put_derivative(colvar, 1, fi) 02705 CALL put_derivative(colvar, 2, fj) 02706 02707 END SUBROUTINE dist_colvar 02708 02709 ! ***************************************************************************** 02713 SUBROUTINE torsion_colvar(colvar,cell,subsys,particles,no_riemann_sheet_op,error) 02714 02715 TYPE(colvar_type), POINTER :: colvar 02716 TYPE(cell_type), POINTER :: cell 02717 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02718 TYPE(particle_type), DIMENSION(:), 02719 OPTIONAL, POINTER :: particles 02720 LOGICAL, INTENT(IN), OPTIONAL :: no_riemann_sheet_op 02721 TYPE(cp_error_type), INTENT(inout) :: error 02722 02723 CHARACTER(len=*), PARAMETER :: routineN = 'torsion_colvar', 02724 routineP = moduleN//':'//routineN 02725 02726 INTEGER :: i, ii 02727 LOGICAL :: failure, no_riemann_sheet 02728 REAL(dp) :: angle, cosine, dedphi, dedxia, dedxib, dedxic, dedxid, dedxt, 02729 dedxu, dedyia, dedyib, dedyic, dedyid, dedyt, dedyu, dedzia, dedzib, 02730 dedzic, dedzid, dedzt, dedzu, dt, e, ftmp(3), o0, rcb, rt2, rtmp(3), 02731 rtru, ru2, sine, ss(3), xba, xca, xcb, xdb, xdc, xt, xtu, xu, yba, yca, 02732 ycb, ydb, ydc, yt, ytu, yu, zba, zca, zcb, zdb, zdc, zt, ztu, zu 02733 REAL(dp), DIMENSION(3, 4) :: rr 02734 TYPE(particle_list_type), POINTER :: particles_i 02735 TYPE(particle_type), DIMENSION(:), 02736 POINTER :: my_particles 02737 02738 failure=.FALSE. 02739 NULLIFY(particles_i) 02740 CPPrecondition(colvar%type_id==torsion_colvar_id,cp_failure_level,routineP,error,failure) 02741 IF (PRESENT(particles)) THEN 02742 my_particles => particles 02743 ELSE 02744 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02745 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02746 my_particles => particles_i%els 02747 END IF 02748 no_riemann_sheet = .FALSE. 02749 IF (PRESENT(no_riemann_sheet_op)) no_riemann_sheet = no_riemann_sheet_op 02750 DO ii=1,4 02751 i=colvar%torsion_param%i_at_tors(ii) 02752 CALL get_coordinates(colvar, i, rtmp, my_particles) 02753 rr(:,ii)=rtmp(1:3) 02754 ENDDO 02755 o0 = colvar%torsion_param%o0 02756 ! ba 02757 ss=MATMUL(cell%h_inv,rr(:,2)-rr(:,1)) 02758 ss=ss-NINT(ss) 02759 ss=MATMUL(cell%hmat,ss) 02760 xba = ss(1) 02761 yba = ss(2) 02762 zba = ss(3) 02763 ! cb 02764 ss=MATMUL(cell%h_inv,rr(:,3)-rr(:,2)) 02765 ss=ss-NINT(ss) 02766 ss=MATMUL(cell%hmat,ss) 02767 xcb = ss(1) 02768 ycb = ss(2) 02769 zcb = ss(3) 02770 ! dc 02771 ss=MATMUL(cell%h_inv,rr(:,4)-rr(:,3)) 02772 ss=ss-NINT(ss) 02773 ss=MATMUL(cell%hmat,ss) 02774 xdc = ss(1) 02775 ydc = ss(2) 02776 zdc = ss(3) 02777 ! 02778 xt = yba*zcb - ycb*zba 02779 yt = zba*xcb - zcb*xba 02780 zt = xba*ycb - xcb*yba 02781 xu = ycb*zdc - ydc*zcb 02782 yu = zcb*xdc - zdc*xcb 02783 zu = xcb*ydc - xdc*ycb 02784 xtu = yt*zu - yu*zt 02785 ytu = zt*xu - zu*xt 02786 ztu = xt*yu - xu*yt 02787 rt2 = xt*xt + yt*yt + zt*zt 02788 ru2 = xu*xu + yu*yu + zu*zu 02789 rtru = SQRT(rt2 * ru2) 02790 IF (rtru .NE. 0.0_dp) THEN 02791 rcb = SQRT(xcb*xcb + ycb*ycb + zcb*zcb) 02792 cosine = (xt*xu + yt*yu + zt*zu) / rtru 02793 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 02794 cosine = MIN(1.0_dp,MAX(-1.0_dp,cosine)) 02795 angle = ACOS(cosine) 02796 IF (sine .LT. 0.0_dp) angle = -angle 02797 ! 02798 dt = angle ! [rad] 02799 dt=MOD(2.0E4_dp*pi+dt-o0,2.0_dp*pi) 02800 IF ( dt .GT. pi ) dt = dt - 2.0_dp*pi 02801 dt=o0+dt 02802 colvar%torsion_param%o0 = dt 02803 ! 02804 ! calculate improper energy and master chain rule term 02805 ! 02806 e = dt 02807 dedphi = 1.0_dp 02808 ! 02809 ! chain rule terms for first derivative components 02810 ! 02811 ! ca 02812 ss=MATMUL(cell%h_inv,rr(:,3)-rr(:,1)) 02813 ss=ss-NINT(ss) 02814 ss=MATMUL(cell%hmat,ss) 02815 xca = ss(1) 02816 yca = ss(2) 02817 zca = ss(3) 02818 ! db 02819 ss=MATMUL(cell%h_inv,rr(:,4)-rr(:,2)) 02820 ss=ss-NINT(ss) 02821 ss=MATMUL(cell%hmat,ss) 02822 xdb = ss(1) 02823 ydb = ss(2) 02824 zdb = ss(3) 02825 ! 02826 dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb) 02827 dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb) 02828 dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb) 02829 dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb) 02830 dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb) 02831 dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb) 02832 ! 02833 ! compute first derivative components for this angle 02834 ! 02835 dedxia = zcb*dedyt - ycb*dedzt 02836 dedyia = xcb*dedzt - zcb*dedxt 02837 dedzia = ycb*dedxt - xcb*dedyt 02838 dedzia = ycb*dedxt - xcb*dedyt 02839 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu 02840 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu 02841 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu 02842 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu 02843 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu 02844 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu 02845 dedxid = zcb*dedyu - ycb*dedzu 02846 dedyid = xcb*dedzu - zcb*dedxu 02847 dedzid = ycb*dedxu - xcb*dedyu 02848 ENDIF 02849 ! 02850 colvar%ss=e 02851 IF (no_riemann_sheet) colvar%ss = ATAN2(SIN(e),COS(e)) 02852 ftmp(1)=dedxia 02853 ftmp(2)=dedyia 02854 ftmp(3)=dedzia 02855 CALL put_derivative(colvar, 1, ftmp) 02856 ftmp(1)=dedxib 02857 ftmp(2)=dedyib 02858 ftmp(3)=dedzib 02859 CALL put_derivative(colvar, 2, ftmp) 02860 ftmp(1)=dedxic 02861 ftmp(2)=dedyic 02862 ftmp(3)=dedzic 02863 CALL put_derivative(colvar, 3, ftmp) 02864 ftmp(1)=dedxid 02865 ftmp(2)=dedyid 02866 ftmp(3)=dedzid 02867 CALL put_derivative(colvar, 4, ftmp) 02868 END SUBROUTINE torsion_colvar 02869 02870 ! ***************************************************************************** 02873 SUBROUTINE qparm_colvar(colvar,cell,subsys,particles,error) 02874 TYPE(colvar_type), POINTER :: colvar 02875 TYPE(cell_type), POINTER :: cell 02876 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 02877 TYPE(particle_type), DIMENSION(:), 02878 OPTIONAL, POINTER :: particles 02879 TYPE(cp_error_type), INTENT(inout) :: error 02880 02881 CHARACTER(len=*), PARAMETER :: routineN = 'qparm_colvar', 02882 routineP = moduleN//':'//routineN 02883 02884 COMPLEX(KIND=dp) :: ylm, ylm_calc, ylm_calc2 02885 INTEGER :: i, ii, j, jj, l, lmm, lpm, m, 02886 n_atoms_from, n_atoms_to 02887 LOGICAL :: failure 02888 REAL(KIND=dp) :: alpha, bond, costheta, dplm, dylm, exp_fac, fact, fi, 02889 ftmp(3), inv_n_atoms_from, nbond, plm, pre_fac, qlm, qlm2, qparm, r12, 02890 r_tmp, rcut, ss(3), theta, x_tmp, xij(3), y_tmp, ymag, ymag1, z_tmp 02891 REAL(KIND=dp), DIMENSION(3) :: grad, grad_nb, grad_nb_tot, 02892 grad_tot, xpi, xpj 02893 TYPE(particle_list_type), POINTER :: particles_i 02894 TYPE(particle_type), DIMENSION(:), 02895 POINTER :: my_particles 02896 02897 failure=.FALSE. 02898 n_atoms_to=colvar%qparm_param%n_atoms_to 02899 n_atoms_from=colvar%qparm_param%n_atoms_from 02900 rcut=colvar%qparm_param%rcut 02901 l=colvar%qparm_param%l 02902 alpha=colvar%qparm_param%alpha 02903 NULLIFY(particles_i) 02904 CPPrecondition(colvar%type_id==qparm_colvar_id,cp_failure_level,routineP,error,failure) 02905 IF (PRESENT(particles)) THEN 02906 my_particles => particles 02907 ELSE 02908 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 02909 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 02910 my_particles => particles_i%els 02911 END IF 02912 02913 qparm=0.0_dp 02914 inv_n_atoms_from=1.0_dp/REAL(n_atoms_from) 02915 DO ii=1,n_atoms_from 02916 i=colvar%qparm_param%i_at_from(ii) 02917 CALL get_coordinates(colvar, i, xpi, my_particles) 02918 qlm = 0.0_dp 02919 qlm2 = 0.0_dp 02920 grad_tot(:) = 0.0_dp 02921 grad_nb_tot(:) = 0.0_dp 02922 DO m = 0, l 02923 ylm = 0.0_dp 02924 grad(:) = 0.0_dp 02925 grad_nb(:) = 0.0_dp 02926 nbond = 0 02927 DO jj=1,n_atoms_to 02928 j=colvar%qparm_param%i_at_to(jj) 02929 IF (i==j) CYCLE 02930 CALL get_coordinates(colvar, j, xpj, my_particles) 02931 ss=MATMUL(cell%h_inv,xpj(:)-xpi(:)) 02932 ss=ss-NINT(ss) 02933 xij=MATMUL(cell%hmat,ss) 02934 r12=SQRT(DOT_PRODUCT(xij,xij)) 02935 IF(r12 < rcut) nbond = nbond + 1.0_dp 02936 x_tmp = xij(1) 02937 y_tmp = xij(2) 02938 z_tmp = xij(3) 02939 r_tmp = r12 02940 IF (ABS(x_tmp) .GT. 1.0E-8_dp) THEN 02941 fi = ATAN(y_tmp/x_tmp) 02942 ELSE 02943 fi = 0.5_dp*pi 02944 ENDIF 02945 costheta = z_tmp/r_tmp 02946 IF(costheta > 1.0_dp) costheta = 1.0_dp 02947 IF(costheta < -1.0_dp) costheta = -1.0_dp 02948 theta = ACOS(costheta) 02949 plm = legendre (costheta, l, m) 02950 dplm = dlegendre (costheta, l, m) 02951 IF ((l+m) > maxfac) THEN 02952 CALL stop_program(routineN,moduleN,__LINE__,& 02953 "(l-m) > maxfac") 02954 END IF 02955 lmm = fac(l-m) 02956 lpm = fac(l+m) 02957 bond = 1.0_dp/(1.0_dp + EXP(alpha*(r_tmp - rcut))) 02958 IF (bond > 1.0_dp) THEN 02959 CALL stop_program(routineN,moduleN,__LINE__,& 02960 "bond > 1.0_dp") 02961 END IF 02962 pre_fac = bond*SQRT(((2*l+1)*lmm)/(4*pi*lpm)) 02963 IF (plm < 0.0_dp) THEN 02964 dylm = -pre_fac*dplm 02965 ELSE 02966 dylm = pre_fac*dplm 02967 ENDIF 02968 ylm_calc = pre_fac*plm*(CMPLX(COS(m*fi),SIN(m*fi))) 02969 ylm = ylm + ylm_calc 02970 ylm_calc2 = pre_fac/bond*plm*(CMPLX(COS(m*fi),SIN(m*fi))) 02971 ymag = (SQRT(ylm_calc2*CONJG(ylm_calc2))) 02972 ! Fi = -dUi/dr, but r_tmp = rj - ri, so have double negative 02973 grad(1) = grad(1) - dylm*(z_tmp*x_tmp)/(r_tmp*r_tmp*r_tmp) 02974 grad(2) = grad(2) - dylm*(z_tmp*y_tmp)/(r_tmp*r_tmp*r_tmp) 02975 grad(3) = grad(3) + dylm*((1.0_dp/r_tmp) - (z_tmp**2/r_tmp**3)) 02976 exp_fac = alpha*EXP(alpha*(r_tmp-rcut))/((1.0_dp+EXP(alpha*(r_tmp-rcut)))**2) 02977 grad_nb(1) = grad_nb(1) - ymag*exp_fac*(x_tmp/r_tmp) 02978 grad_nb(2) = grad_nb(2) - ymag*exp_fac*(y_tmp/r_tmp) 02979 grad_nb(3) = grad_nb(3) - ymag*exp_fac*(z_tmp/r_tmp) 02980 ENDDO 02981 IF (m .GT. 0) THEN 02982 fact = 2.0_dp 02983 ELSE 02984 fact = 1.0_dp 02985 ENDIF 02986 IF (nbond < 0.0001_dp) THEN 02987 CALL stop_program(routineN,moduleN,__LINE__,& 02988 "NBOND = 0.0") 02989 END IF 02990 ylm = ylm/nbond 02991 grad(:) = grad(:)/nbond 02992 grad_nb(:) = grad_nb(:)/nbond 02993 ymag1 = ylm*CONJG(ylm) 02994 qlm = qlm + fact*ymag1 02995 grad_tot(:) = grad_tot(:) + fact*grad(:) 02996 grad_nb_tot(:) = grad_nb_tot(:) + fact*grad_nb(:) 02997 ENDDO 02998 pre_fac = (4.0_dp*pi)/(2.0_dp*l + 1) 02999 qlm = SQRT(pre_fac*qlm) 03000 qlm2 = SQRT(pre_fac*qlm2) 03001 qparm=qparm+qlm 03002 ftmp(:) = SQRT(pre_fac)*(grad_tot(:) + grad_nb_tot(:)) 03003 CALL put_derivative(colvar, ii, ftmp) 03004 ENDDO 03005 colvar%ss=qparm*inv_n_atoms_from 03006 colvar%dsdr(:,:) = colvar%dsdr(:,:)*inv_n_atoms_from 03007 03008 END SUBROUTINE qparm_colvar 03009 03010 ! ***************************************************************************** 03016 SUBROUTINE hydronium_colvar(colvar,cell,subsys,particles,error) 03017 TYPE(colvar_type), POINTER :: colvar 03018 TYPE(cell_type), POINTER :: cell 03019 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 03020 TYPE(particle_type), DIMENSION(:), 03021 OPTIONAL, POINTER :: particles 03022 TYPE(cp_error_type), INTENT(inout), 03023 OPTIONAL :: error 03024 03025 CHARACTER(len=*), PARAMETER :: routineN = 'hydronium_colvar', 03026 routineP = moduleN//':'//routineN 03027 03028 INTEGER :: i, ii, jj, n_hydrogens, 03029 n_oxygens, p, pnh, pno, q, 03030 qnh, qno, stat 03031 LOGICAL :: failure 03032 REAL(dp) :: fscalar, invden, lambda, nh, 03033 num, qtot, r12, r_OH, r_OO, 03034 rdist, ss(3), xij(3) 03035 REAL(dp), DIMENSION(3) :: xpi, xpj 03036 REAL(dp), DIMENSION(:) :: M, nhcoord, no, qloc 03037 03038 ALLOCATABLE :: nhcoord,M,& 03039 qloc,no 03040 REAL(dp), DIMENSION(:,:,:) :: dfunc_nh,dfunc_no,dM 03041 ALLOCATABLE :: dfunc_nh,dfunc_no,dM 03042 TYPE(particle_list_type), POINTER :: particles_i 03043 TYPE(particle_type), DIMENSION(:), 03044 POINTER :: my_particles 03045 03046 failure=.FALSE. 03047 n_oxygens=colvar%hydronium_param%n_oxygens 03048 n_hydrogens=colvar%hydronium_param%n_hydrogens 03049 nh=colvar%hydronium_param%nh 03050 pnh=colvar%hydronium_param%pnh 03051 qnh=colvar%hydronium_param%qnh 03052 pno=colvar%hydronium_param%pno 03053 qno=colvar%hydronium_param%qno 03054 r_OO=colvar%hydronium_param%r_OO 03055 r_OH=colvar%hydronium_param%r_OH 03056 lambda=colvar%hydronium_param%lambda 03057 p=colvar%hydronium_param%p 03058 q=colvar%hydronium_param%q 03059 03060 NULLIFY(particles_i) 03061 CPPrecondition(colvar%type_id==hydronium_colvar_id,cp_failure_level,routineP,error,failure) 03062 IF (PRESENT(particles)) THEN 03063 my_particles => particles 03064 ELSE 03065 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 03066 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 03067 my_particles => particles_i%els 03068 END IF 03069 03070 ALLOCATE ( dfunc_nh ( 3, n_oxygens, n_hydrogens ), stat=stat ) 03071 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03072 ALLOCATE ( nhcoord ( n_oxygens), stat=stat ) 03073 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03074 ALLOCATE ( M ( n_oxygens ), stat=stat ) 03075 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03076 ALLOCATE ( dM ( 3, n_oxygens,n_hydrogens ), stat=stat ) 03077 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03078 03079 ALLOCATE ( dfunc_no ( 3, n_oxygens, n_oxygens ), stat=stat ) 03080 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03081 ALLOCATE ( no ( n_oxygens), stat=stat ) 03082 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03083 03084 ALLOCATE ( qloc ( n_oxygens), stat=stat ) 03085 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03086 03087 ! Zero Arrays: 03088 dfunc_nh = 0._dp 03089 dfunc_no = 0._dp 03090 M = 0._dp 03091 dM = 0._dp 03092 no = 0._dp 03093 qloc = 0._dp 03094 nhcoord = 0._dp 03095 DO ii=1,n_oxygens 03096 i=colvar%hydronium_param%i_oxygens(ii) 03097 xpi(:)=my_particles(i)%r(1:3) 03098 ! Computing M( n ( ii ) ) 03099 DO jj=1,n_hydrogens 03100 i=colvar%hydronium_param%i_hydrogens(jj) 03101 xpj(:)=my_particles(i)%r(1:3) 03102 ss=MATMUL(cell%h_inv,xpi(:)-xpj(:)) 03103 ss=ss-NINT(ss) 03104 xij=MATMUL(cell%hmat,ss) 03105 r12=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 03106 rdist = r12/r_OH 03107 num=(1.0_dp-rdist**pnh) 03108 invden=1.0_dp/(1.0_dp-rdist**qnh) 03109 IF ( ABS ( invden ) < 1.e-10_dp ) invden = 1.e-10_dp 03110 fscalar = ((-pnh*(rdist**(pnh-1))*invden) & 03111 + num*(invden)**2*qnh*(rdist**(qnh-1)))/(r12*r_OH) 03112 03113 dfunc_nh(1,ii,jj)= xij (1) * fscalar 03114 dfunc_nh(2,ii,jj)= xij (2) * fscalar 03115 dfunc_nh(3,ii,jj)= xij (3) * fscalar 03116 03117 nhcoord(ii)=nhcoord(ii) + num*invden 03118 END DO 03119 M(ii) = 1.0_dp-(1.0_dp - (nhcoord(ii)/nh)**p)/& 03120 (1.0_dp - (nhcoord(ii)/nh)**q) 03121 03122 ! Computing no ( ii ) 03123 DO jj=1,n_oxygens 03124 i=colvar%hydronium_param%i_oxygens(jj) 03125 xpj(:)=my_particles(i)%r(1:3) 03126 ss=MATMUL(cell%h_inv,xpi(:)-xpj(:)) 03127 ss=ss-NINT(ss) 03128 xij=MATMUL(cell%hmat,ss) 03129 r12=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 03130 IF ( r12 < 1.e-3_dp ) CYCLE 03131 rdist = r12/r_OO 03132 num=(1.0_dp-rdist**pno) 03133 invden=1.0_dp/(1.0_dp-rdist**qno) 03134 IF ( ABS ( invden ) < 1.e-10_dp ) invden = 1.e-10_dp 03135 fscalar = ((-pno*(rdist**(pno-1))*invden) & 03136 + num*(invden)**2*qno*(rdist**(qno-1)))/(r12*r_OO) 03137 03138 dfunc_no(1,ii,jj)= xij(1)*fscalar 03139 dfunc_no(2,ii,jj)= xij(2)*fscalar 03140 dfunc_no(3,ii,jj)= xij(3)*fscalar 03141 03142 no(ii)=no(ii)+ num*invden 03143 END DO 03144 END DO 03145 03146 ! computing qloc and Q 03147 qtot = 0._dp 03148 DO ii = 1, n_oxygens 03149 qloc ( ii ) = EXP ( lambda * M (ii) * no ( ii ) ) 03150 qtot = qtot + qloc ( ii ) 03151 END DO 03152 ! compute forces 03153 DO ii = 1, n_oxygens 03154 ! Computing f_OH 03155 DO jj=1,n_hydrogens 03156 dM(1,ii,jj)=(p*((nhcoord(ii)/nh)**(p-1))*dfunc_nh(1,ii,jj))/nh/ & 03157 (1.0_dp - (nhcoord(ii)/nh)**q) - & 03158 (1.0_dp - (nhcoord(ii)/nh)**p)/ & 03159 ((1.0_dp - (nhcoord(ii)/nh)**q)**2) *& 03160 q*dfunc_nh(1,ii,jj)*(nhcoord(ii)/nh)**(q-1)/nh 03161 dM(2,ii,jj)=(p*((nhcoord(ii)/nh)**(p-1))*dfunc_nh(2,ii,jj))/nh/ & 03162 (1.0_dp - (nhcoord(ii)/nh)**q) - & 03163 (1.0_dp - (nhcoord(ii)/nh)**p)/ & 03164 ((1.0_dp - (nhcoord(ii)/nh)**q)**2) *& 03165 q*dfunc_nh(2,ii,jj)*(nhcoord(ii)/nh)**(q-1)/nh 03166 dM(3,ii,jj)=(p*((nhcoord(ii)/nh)**(p-1))*dfunc_nh(3,ii,jj))/nh/ & 03167 (1.0_dp - (nhcoord(ii)/nh)**q) - & 03168 (1.0_dp - (nhcoord(ii)/nh)**p)/ & 03169 ((1.0_dp - (nhcoord(ii)/nh)**q)**2) *& 03170 q*dfunc_nh(3,ii,jj)*(nhcoord(ii)/nh)**(q-1)/nh 03171 03172 colvar%dsdr(1,ii)=colvar%dsdr(1,ii)+qloc(ii)*dM(1,ii,jj)*no(ii)/qtot 03173 colvar%dsdr(1,n_oxygens+jj)=colvar%dsdr(1,n_oxygens+jj) & 03174 -qloc(ii)*dM(1,ii,jj)*no(ii)/qtot 03175 colvar%dsdr(2,ii)=colvar%dsdr(2,ii)+qloc(ii)*dM(2,ii,jj)*no(ii)/qtot 03176 colvar%dsdr(2,n_oxygens+jj)=colvar%dsdr(2,n_oxygens+jj) & 03177 -qloc(ii)*dM(2,ii,jj)*no(ii)/qtot 03178 colvar%dsdr(3,ii)=colvar%dsdr(3,ii)+qloc(ii)*dM(3,ii,jj)*no(ii)/qtot 03179 colvar%dsdr(3,n_oxygens+jj)=colvar%dsdr(3,n_oxygens+jj) & 03180 -qloc(ii)*dM(3,ii,jj)*no(ii)/qtot 03181 END DO 03182 ! Computing f_OO 03183 DO jj=1,n_oxygens 03184 colvar%dsdr(1,ii)=colvar%dsdr(1,ii)+qloc(ii)*M(ii)*dfunc_no(1,ii,jj)/qtot 03185 colvar%dsdr(1,jj)=colvar%dsdr(1,jj) & 03186 -qloc(ii)*M(ii)*dfunc_no(1,ii,jj)/qtot 03187 colvar%dsdr(2,ii)=colvar%dsdr(2,ii)+qloc(ii)*M(ii)*dfunc_no(2,ii,jj)/qtot 03188 colvar%dsdr(2,jj)=colvar%dsdr(2,jj) & 03189 -qloc(ii)*M(ii)*dfunc_no(2,ii,jj)/qtot 03190 colvar%dsdr(3,ii)=colvar%dsdr(3,ii)+qloc(ii)*M(ii)*dfunc_no(3,ii,jj)/qtot 03191 colvar%dsdr(3,jj)=colvar%dsdr(3,jj) & 03192 -qloc(ii)*M(ii)*dfunc_no(3,ii,jj)/qtot 03193 END DO 03194 END DO 03195 03196 colvar%ss=LOG(qtot)/lambda 03197 DEALLOCATE ( dfunc_nh, stat=stat ) 03198 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03199 DEALLOCATE ( nhcoord, stat=stat ) 03200 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03201 DEALLOCATE ( M, stat=stat ) 03202 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03203 DEALLOCATE ( dM, stat=stat ) 03204 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03205 DEALLOCATE ( dfunc_no, stat=stat ) 03206 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03207 DEALLOCATE ( no, stat=stat ) 03208 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03209 DEALLOCATE ( qloc, stat=stat ) 03210 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03211 END SUBROUTINE hydronium_colvar 03212 03213 ! ***************************************************************************** 03219 SUBROUTINE coord_colvar(colvar,cell,subsys,particles,error) 03220 TYPE(colvar_type), POINTER :: colvar 03221 TYPE(cell_type), POINTER :: cell 03222 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 03223 TYPE(particle_type), DIMENSION(:), 03224 OPTIONAL, POINTER :: particles 03225 TYPE(cp_error_type), INTENT(inout) :: error 03226 03227 CHARACTER(len=*), PARAMETER :: routineN = 'coord_colvar', 03228 routineP = moduleN//':'//routineN 03229 03230 INTEGER :: i, ii, j, jj, k, kk, 03231 n_atoms_from, n_atoms_to_a, 03232 n_atoms_to_b, p_a, p_b, q_a, 03233 q_b 03234 LOGICAL :: failure 03235 REAL(dp) :: dfunc_ij, dfunc_jk, func_ij, func_jk, func_k, 03236 inv_n_atoms_from, invden_ij, invden_jk, ncoord, num_ij, num_jk, r_0_a, 03237 r_0_b, rdist_ij, rdist_jk, rij, rjk 03238 REAL(dp), DIMENSION(3) :: ftmp_i, ftmp_j, ftmp_k, ss, 03239 xij, xjk, xpi, xpj, xpk 03240 TYPE(particle_list_type), POINTER :: particles_i 03241 TYPE(particle_type), DIMENSION(:), 03242 POINTER :: my_particles 03243 03244 failure=.FALSE. 03245 ! If we defined the coordination number with KINDS then we have still 03246 ! to fill few missing informations... 03247 NULLIFY(particles_i) 03248 CPPrecondition(colvar%type_id==coord_colvar_id,cp_failure_level,routineP,error,failure) 03249 IF (PRESENT(particles)) THEN 03250 my_particles => particles 03251 ELSE 03252 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 03253 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 03254 my_particles => particles_i%els 03255 END IF 03256 n_atoms_to_a=colvar%coord_param%n_atoms_to 03257 n_atoms_to_b=colvar%coord_param%n_atoms_to_b 03258 n_atoms_from=colvar%coord_param%n_atoms_from 03259 p_a=colvar%coord_param%nncrd 03260 q_a=colvar%coord_param%ndcrd 03261 r_0_a=colvar%coord_param%r_0 03262 p_b=colvar%coord_param%nncrd_b 03263 q_b=colvar%coord_param%ndcrd_b 03264 r_0_b=colvar%coord_param%r_0_b 03265 03266 ncoord=0.0_dp 03267 inv_n_atoms_from=1.0_dp/REAL(n_atoms_from) 03268 DO ii=1,n_atoms_from 03269 i=colvar%coord_param%i_at_from(ii) 03270 CALL get_coordinates(colvar, i, xpi, my_particles) 03271 DO jj=1,n_atoms_to_a 03272 j=colvar%coord_param%i_at_to(jj) 03273 CALL get_coordinates(colvar, j, xpj, my_particles) 03274 ! define coordination of atom A with itself to be 0. also fixes rij==0 for the force calculation 03275 IF (i.EQ.j) CYCLE 03276 ss=MATMUL(cell%h_inv,xpi(:)-xpj(:)) 03277 ss=ss-NINT(ss) 03278 xij=MATMUL(cell%hmat,ss) 03279 rij=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 03280 IF(rij < 1.0e-8_dp)CYCLE 03281 rdist_ij = rij/r_0_a 03282 IF(ABS(1.0_dp-rdist_ij) > EPSILON(0.0_dp)*1.0E+4_dp) THEN 03283 num_ij=(1.0_dp-rdist_ij**p_a) 03284 invden_ij=1.0_dp/(1.0_dp-rdist_ij**q_a) 03285 func_ij=num_ij*invden_ij 03286 IF (rij < 1.0E-8_dp) THEN 03287 ! provide the correct limit of the derivative 03288 dfunc_ij = 0.0_dp 03289 ELSE 03290 dfunc_ij= (- p_a*rdist_ij **(p_a-1) *invden_ij & 03291 + num_ij*(invden_ij)**2 * q_a *rdist_ij **(q_a-1))/(rij*r_0_a) 03292 END IF 03293 ELSE 03294 ! Provide the correct limit for function value and derivative 03295 func_ij = REAL(p_a,KIND=dp)/REAL(q_a,KIND=dp) 03296 dfunc_ij = REAL(p_a,KIND=dp)*REAL((-q_a+p_a),KIND=dp)/(REAL(2*q_a,KIND=dp)*r_0_a) 03297 END IF 03298 IF(n_atoms_to_b /=0) THEN 03299 func_k = 0.0_dp 03300 DO kk = 1,n_atoms_to_b 03301 k=colvar%coord_param%i_at_to_b(kk) 03302 IF (k.EQ.j) CYCLE 03303 CALL get_coordinates(colvar, k, xpk, my_particles) 03304 ss=MATMUL(cell%h_inv,xpj(:)-xpk(:)) 03305 ss=ss-NINT(ss) 03306 xjk=MATMUL(cell%hmat,ss) 03307 rjk=SQRT(xjk(1)**2+xjk(2)**2+xjk(3)**2) 03308 IF(rjk < 1.0e-8_dp)CYCLE 03309 rdist_jk = rjk/r_0_b 03310 IF(ABS(1.0_dp-rdist_jk) > EPSILON(0.0_dp)*1.0E+4_dp) THEN 03311 num_jk=(1.0_dp-rdist_jk**p_b) 03312 invden_jk=1.0_dp/(1.0_dp-rdist_jk**q_b) 03313 func_jk=num_jk*invden_jk 03314 IF (rjk < 1.0E-8_dp) THEN 03315 ! provide the correct limit of the derivative 03316 dfunc_jk = 0.0_dp 03317 ELSE 03318 dfunc_jk= (- p_b*rdist_jk **(p_b-1) *invden_jk & 03319 + num_jk*(invden_jk)**2 * q_b *rdist_jk **(q_b-1))/(rjk*r_0_b) 03320 END IF 03321 ELSE 03322 ! Provide the correct limit for function value and derivative 03323 func_jk = REAL(p_b,KIND=dp)/REAL(q_b,KIND=dp) 03324 dfunc_jk = REAL(p_b,KIND=dp)*REAL((-q_b+p_b),KIND=dp)/(REAL(2*q_b,KIND=dp)*r_0_b) 03325 ENDIF 03326 func_k = func_k + func_jk 03327 ftmp_k = - func_ij * dfunc_jk * xjk 03328 CALL put_derivative(colvar, n_atoms_from+n_atoms_to_a+kk, ftmp_k) 03329 03330 ftmp_j = - dfunc_ij * xij * func_jk + func_ij * dfunc_jk * xjk 03331 CALL put_derivative(colvar, n_atoms_from+jj, ftmp_j) 03332 END DO 03333 ELSE 03334 func_k = 1.0_dp 03335 dfunc_jk = 0.0_dp 03336 ftmp_j = -dfunc_ij*xij 03337 CALL put_derivative(colvar, n_atoms_from+jj, ftmp_j) 03338 END IF 03339 ncoord=ncoord+func_ij*func_k 03340 ftmp_i = dfunc_ij*xij*func_k 03341 CALL put_derivative(colvar, ii, ftmp_i) 03342 ENDDO 03343 ENDDO 03344 colvar%ss=ncoord*inv_n_atoms_from 03345 colvar%dsdr(:,:) = colvar%dsdr(:,:)*inv_n_atoms_from 03346 END SUBROUTINE coord_colvar 03347 03348 SUBROUTINE mindist_colvar(colvar,cell, subsys, particles, error) 03349 03350 TYPE(colvar_type), POINTER :: colvar 03351 TYPE(cell_type), POINTER :: cell 03352 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 03353 TYPE(particle_type), DIMENSION(:), 03354 OPTIONAL, POINTER :: particles 03355 TYPE(cp_error_type), INTENT(inout) :: error 03356 03357 CHARACTER(len=*), PARAMETER :: routineN = 'mindist_colvar', 03358 routineP = moduleN//':'//routineN 03359 03360 INTEGER :: i, ii, j, jj, n_coord_from, 03361 n_coord_to, n_dist_from, p, 03362 q, stat 03363 LOGICAL :: failure 03364 REAL(dp) :: den_n, den_Q, fscalar, ftemp_i(3), inv_den_n, inv_den_Q, 03365 lambda, num_n, num_Q, Qfunc, r12, r_cut, rfact, rij(3), rpi(3), rpj(3) 03366 REAL(dp), DIMENSION(:), POINTER :: dqfunc_dnL, expnL, nLcoord, 03367 sum_rij 03368 REAL(dp), DIMENSION(:, :, :), POINTER :: dnLcoord, dqfunc_dr 03369 TYPE(particle_list_type), POINTER :: particles_i 03370 TYPE(particle_type), DIMENSION(:), 03371 POINTER :: my_particles 03372 03373 failure=.FALSE. 03374 ! If we defined the coordination number with KINDS then we have still 03375 ! to fill few missing informations... 03376 NULLIFY(particles_i) 03377 CPPrecondition(colvar%type_id==mindist_colvar_id,cp_failure_level,routineP,error,failure) 03378 IF (PRESENT(particles)) THEN 03379 my_particles => particles 03380 ELSE 03381 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 03382 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 03383 my_particles => particles_i%els 03384 END IF 03385 03386 n_dist_from=colvar%mindist_param%n_dist_from 03387 n_coord_from=colvar%mindist_param%n_coord_from 03388 n_coord_to=colvar%mindist_param%n_coord_to 03389 p=colvar%mindist_param%p_exp 03390 q=colvar%mindist_param%q_exp 03391 r_cut=colvar%mindist_param%r_cut 03392 lambda=colvar%mindist_param%lambda 03393 03394 03395 NULLIFY(nLcoord,dnLcoord,dqfunc_dr,dqfunc_dnL,expnL,sum_rij) 03396 ALLOCATE (nLcoord(n_coord_from),stat=stat) 03397 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03398 ALLOCATE (dnLcoord(3,n_coord_from,n_coord_to),stat=stat) 03399 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03400 ALLOCATE (expnL(n_coord_from),stat=stat) 03401 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03402 ALLOCATE (sum_rij(n_coord_from),stat=stat) 03403 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03404 ALLOCATE (dqfunc_dr(3,n_dist_from,n_coord_from),stat=stat) 03405 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03406 ALLOCATE (dqfunc_dnL(n_coord_from),stat=stat) 03407 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03408 03409 ! coordination numbers 03410 nLcoord = 0.0_dp 03411 dnLcoord = 0.0_dp 03412 expnL = 0.0_dp 03413 den_Q = 0.0_dp 03414 DO i = 1,n_coord_from 03415 ii = colvar%mindist_param%i_coord_from(i) 03416 rpi=my_particles(ii)%r(1:3) 03417 DO j = 1,n_coord_to 03418 jj = colvar%mindist_param%i_coord_to(j) 03419 rpj=my_particles(jj)%r(1:3) 03420 rij = pbc(rpj,rpi,cell) 03421 r12 = SQRT(rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)) 03422 rfact=r12/r_cut 03423 num_n = 1.0_dp-rfact**p 03424 den_n = 1.0_dp-rfact**q 03425 inv_den_n = 1.0_dp/den_n 03426 IF(ABS(inv_den_n)<1.e-10_dp ) THEN 03427 inv_den_n=1.e-10_dp 03428 num_n = ABS(num_n) 03429 END IF 03430 03431 fscalar = (-p*rfact**(p-1)+num_n*q*rfact**(q-1)*inv_den_n)*inv_den_n/(r_cut*r12) 03432 03433 dnLcoord(1,i,j) = rij(1)*fscalar 03434 dnLcoord(2,i,j) = rij(2)*fscalar 03435 dnLcoord(3,i,j) = rij(3)*fscalar 03436 03437 nLcoord(i) = nLcoord(i) + num_n*inv_den_n 03438 END DO 03439 expnL(i) = EXP(lambda*nLcoord(i)) 03440 !dbg 03441 ! write(*,*) ii,nLcoord(i),expnL(i) 03442 !dbg 03443 den_Q = den_Q + expnL(i) 03444 END DO 03445 inv_den_Q = 1.0_dp/den_Q 03446 03447 qfunc = 0.0_dp 03448 dqfunc_dr = 0.0_dp 03449 dqfunc_dnL = 0.0_dp 03450 num_Q = 0.0_dp 03451 sum_rij = 0.0_dp 03452 DO i = 1,n_dist_from 03453 ii = colvar%mindist_param%i_dist_from(i) 03454 rpi=my_particles(ii)%r(1:3) 03455 DO j = 1,n_coord_from 03456 jj = colvar%mindist_param%i_coord_from(j) 03457 rpj=my_particles(jj)%r(1:3) 03458 rij = pbc(rpj,rpi,cell) 03459 r12 = SQRT(rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)) 03460 03461 !dbg 03462 ! write(*,*) ii,jj,rpi(1:3),rpj(1:3),rij(1:3),r12 03463 !dbg 03464 num_Q = num_Q + r12 * expnL(j) 03465 03466 sum_rij(j) = sum_rij(j) + r12 03467 dqfunc_dr(1,i,j) = expnL(j)* rij(1)/r12 03468 dqfunc_dr(2,i,j) = expnL(j)* rij(2)/r12 03469 dqfunc_dr(3,i,j) = expnL(j)* rij(3)/r12 03470 03471 END DO 03472 03473 END DO 03474 03475 ! Function and derivatives 03476 qfunc = num_Q * inv_den_Q 03477 dqfunc_dr = dqfunc_dr * inv_den_Q 03478 colvar%ss=qfunc 03479 !dbg 03480 ! write(*,*) ' ss ', colvar%ss 03481 ! stop 03482 !dbg 03483 03484 DO i = 1,n_coord_from 03485 dqfunc_dnL(i) = lambda * expnL(i) * inv_den_Q*( sum_rij(i) - num_Q * inv_den_Q) 03486 END DO 03487 03488 !Compute Forces 03489 DO i = 1,n_dist_from 03490 DO j = 1,n_coord_from 03491 ftemp_i(1) = dqfunc_dr(1,i,j) 03492 ftemp_i(2) = dqfunc_dr(2,i,j) 03493 ftemp_i(3) = dqfunc_dr(3,i,j) 03494 03495 CALL put_derivative(colvar,i,ftemp_i) 03496 CALL put_derivative(colvar,j+n_dist_from,-ftemp_i) 03497 03498 END DO 03499 END DO 03500 DO i = 1,n_coord_from 03501 DO j = 1,n_coord_to 03502 ftemp_i(1) = dqfunc_dnL(i)*dnLcoord(1,i,j) 03503 ftemp_i(2) = dqfunc_dnL(i)*dnLcoord(2,i,j) 03504 ftemp_i(3) = dqfunc_dnL(i)*dnLcoord(3,i,j) 03505 03506 CALL put_derivative(colvar,i+n_dist_from,ftemp_i) 03507 CALL put_derivative(colvar,j+n_dist_from+n_coord_from,-ftemp_i) 03508 03509 END DO 03510 END DO 03511 03512 03513 03514 03515 DEALLOCATE(nLcoord,stat=stat) 03516 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03517 DEALLOCATE(dnLcoord,stat=stat) 03518 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03519 DEALLOCATE (expnL,stat=stat) 03520 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03521 DEALLOCATE (dqfunc_dr,stat=stat) 03522 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03523 DEALLOCATE (sum_rij,stat=stat) 03524 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03525 DEALLOCATE (dqfunc_dnL,stat=stat) 03526 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03527 03528 END SUBROUTINE mindist_colvar 03529 03530 ! ***************************************************************************** 03534 SUBROUTINE combine_colvar(colvar,cell,subsys,particles,error) 03535 TYPE(colvar_type), POINTER :: colvar 03536 TYPE(cell_type), POINTER :: cell 03537 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 03538 TYPE(particle_type), DIMENSION(:), 03539 OPTIONAL, POINTER :: particles 03540 TYPE(cp_error_type), INTENT(inout) :: error 03541 03542 CHARACTER(len=*), PARAMETER :: routineN = 'combine_colvar', 03543 routineP = moduleN//':'//routineN 03544 03545 CHARACTER(LEN=default_string_length) :: def_error, this_error 03546 CHARACTER(LEN=default_string_length), 03547 ALLOCATABLE, DIMENSION(:) :: my_par 03548 INTEGER :: i, ii, j, ncolv, ndim, stat 03549 LOGICAL :: failure 03550 REAL(dp) :: err 03551 REAL(dp), ALLOCATABLE, DIMENSION(:) :: dss_vals, my_val, ss_vals 03552 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fi 03553 TYPE(particle_list_type), POINTER :: particles_i 03554 TYPE(particle_type), DIMENSION(:), 03555 POINTER :: my_particles 03556 03557 failure=.FALSE. 03558 CPPrecondition(colvar%type_id==combine_colvar_id,cp_failure_level,routineP,error,failure) 03559 IF (PRESENT(particles)) THEN 03560 my_particles => particles 03561 ELSE 03562 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 03563 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 03564 my_particles => particles_i%els 03565 END IF 03566 03567 ncolv=SIZE(colvar%combine_cvs_param%colvar_p) 03568 ALLOCATE(ss_vals(ncolv),stat=stat) 03569 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03570 ALLOCATE(dss_vals(ncolv),stat=stat) 03571 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03572 03573 ! Evaluate the individual COLVARs 03574 DO i=1,ncolv 03575 CALL colvar_recursive_eval(colvar%combine_cvs_param%colvar_p(i)%colvar,cell,my_particles,error) 03576 ss_vals(i)=colvar%combine_cvs_param%colvar_p(i)%colvar%ss 03577 ENDDO 03578 03579 ! Evaluate the combination of the COLVARs 03580 CALL initf(1) 03581 ndim = SIZE(colvar%combine_cvs_param%c_parameters)+& 03582 SIZE(colvar%combine_cvs_param%variables) 03583 ALLOCATE(my_par(ndim),stat=stat) 03584 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03585 my_par(1:SIZE(colvar%combine_cvs_param%variables)) =colvar%combine_cvs_param%variables 03586 my_par(SIZE(colvar%combine_cvs_param%variables)+1:)=colvar%combine_cvs_param%c_parameters 03587 ALLOCATE(my_val(ndim),stat=stat) 03588 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03589 my_val(1:SIZE(colvar%combine_cvs_param%variables)) =ss_vals 03590 my_val(SIZE(colvar%combine_cvs_param%variables)+1:)=colvar%combine_cvs_param%v_parameters 03591 CALL parsef(1,TRIM(colvar%combine_cvs_param%function),my_par) 03592 colvar%ss=evalf(1,my_val) 03593 DO i = 1, ncolv 03594 dss_vals(i) = evalfd(1,i,my_val,colvar%combine_cvs_param%dx, err) 03595 IF ((ABS(err)>colvar%combine_cvs_param%lerr)) THEN 03596 WRITE(this_error,"(A,G12.6,A)")"(",err,")" 03597 WRITE(def_error,"(A,G12.6,A)")"(",colvar%combine_cvs_param%lerr,")" 03598 CALL compress(this_error,.TRUE.) 03599 CALL compress(def_error,.TRUE.) 03600 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 03601 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)//& 03602 ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'//& 03603 TRIM(def_error)//' . '//& 03604 CPSourceFileRef,& 03605 only_ionode=.TRUE.) 03606 END IF 03607 END DO 03608 DEALLOCATE(my_val,stat=stat) 03609 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03610 DEALLOCATE(my_par,stat=stat) 03611 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03612 CALL finalizef() 03613 03614 ! Evaluate forces 03615 ALLOCATE(fi(3,colvar%n_atom_s),stat=stat) 03616 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03617 ii=0 03618 DO i=1,ncolv 03619 DO j=1,colvar%combine_cvs_param%colvar_p(i)%colvar%n_atom_s 03620 ii=ii+1 03621 fi(:,ii)=colvar%combine_cvs_param%colvar_p(i)%colvar%dsdr(:,j)*dss_vals(i) 03622 END DO 03623 END DO 03624 03625 DO i=1,colvar%n_atom_s 03626 CALL put_derivative(colvar,i,fi(:,i)) 03627 END DO 03628 03629 DEALLOCATE(fi,stat=stat) 03630 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03631 DEALLOCATE(ss_vals,stat=stat) 03632 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03633 DEALLOCATE(dss_vals,stat=stat) 03634 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03635 END SUBROUTINE combine_colvar 03636 03637 ! ***************************************************************************** 03647 SUBROUTINE reaction_path_colvar(colvar,cell,subsys,particles,error) 03648 TYPE(colvar_type), POINTER :: colvar 03649 TYPE(cell_type), POINTER :: cell 03650 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 03651 TYPE(particle_type), DIMENSION(:), 03652 OPTIONAL, POINTER :: particles 03653 TYPE(cp_error_type), INTENT(inout) :: error 03654 03655 CHARACTER(len=*), PARAMETER :: routineN = 'reaction_path_colvar', 03656 routineP = moduleN//':'//routineN 03657 03658 LOGICAL :: failure 03659 TYPE(particle_list_type), POINTER :: particles_i 03660 TYPE(particle_type), DIMENSION(:), 03661 POINTER :: my_particles 03662 03663 failure=.FALSE. 03664 03665 CPPrecondition(colvar%type_id==reaction_path_colvar_id,cp_failure_level,routineP,error,failure) 03666 IF (PRESENT(particles)) THEN 03667 my_particles => particles 03668 ELSE 03669 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 03670 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 03671 my_particles => particles_i%els 03672 END IF 03673 03674 IF(colvar%reaction_path_param%dist_rmsd) THEN 03675 CALL rpath_dist_rmsd(colvar,my_particles,error=error) 03676 ELSEIF(colvar%reaction_path_param%rmsd) THEN 03677 CALL rpath_rmsd(colvar,my_particles,error=error) 03678 ELSE 03679 CALL rpath_colvar(colvar,cell,my_particles,error) 03680 END IF 03681 03682 END SUBROUTINE reaction_path_colvar 03683 03684 ! ***************************************************************************** 03690 SUBROUTINE rpath_colvar(colvar,cell,particles,error) 03691 TYPE(colvar_type), POINTER :: colvar 03692 TYPE(cell_type), POINTER :: cell 03693 TYPE(particle_type), DIMENSION(:), 03694 POINTER :: particles 03695 TYPE(cp_error_type), INTENT(inout) :: error 03696 03697 CHARACTER(len=*), PARAMETER :: routineN = 'rpath_colvar', 03698 routineP = moduleN//':'//routineN 03699 03700 INTEGER :: i, iend, ii, istart, j, k, 03701 ncolv, nconf, stat 03702 LOGICAL :: failure 03703 REAL(dp) :: lambda, step_size 03704 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s1, ss_vals 03705 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ds1, f_vals, fi, s1v 03706 REAL(dp), ALLOCATABLE, 03707 DIMENSION(:, :, :) :: ds1v 03708 03709 istart=colvar%reaction_path_param%function_bounds(1) 03710 iend=colvar%reaction_path_param%function_bounds(2) 03711 03712 nconf=colvar%reaction_path_param%nr_frames 03713 step_size=colvar%reaction_path_param%step_size 03714 ncolv=colvar%reaction_path_param%n_components 03715 lambda=colvar%reaction_path_param%lambda 03716 ALLOCATE(f_vals(ncolv,istart:iend),stat=stat) 03717 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03718 f_vals=colvar%reaction_path_param%f_vals 03719 ALLOCATE(ss_vals(ncolv),stat=stat) 03720 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03721 03722 DO i=1,ncolv 03723 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar,cell,particles,error) 03724 ss_vals(i)=colvar%reaction_path_param%colvar_p(i)%colvar%ss 03725 ENDDO 03726 03727 ALLOCATE(s1v(2,istart:iend),stat=stat) 03728 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03729 ALLOCATE(ds1v(ncolv,2,istart:iend),stat=stat) 03730 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03731 03732 ALLOCATE(s1(2),stat=stat) 03733 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03734 ALLOCATE(ds1(ncolv,2),stat=stat) 03735 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03736 03737 DO k=istart,iend 03738 s1v(1,k)=REAL(k,kind=dp)*step_size*EXP(-lambda*DOT_PRODUCT(ss_vals(:)-f_vals(:,k),ss_vals(:)-f_vals(:,k))) 03739 s1v(2,k)=EXP(-lambda*DOT_PRODUCT(ss_vals(:)-f_vals(:,k),ss_vals(:)-f_vals(:,k))) 03740 DO j=1,ncolv 03741 ds1v(j,1,k) = f_vals(j,k)* s1v(1,k) 03742 ds1v(j,2,k) = f_vals(j,k) * s1v(2,k) 03743 END DO 03744 END DO 03745 DO i=1,2 03746 s1(i)=accurate_sum(s1v(i,:)) 03747 DO j=1,ncolv 03748 ds1(j,i)=accurate_sum(ds1v(j,i,:)) 03749 END DO 03750 END DO 03751 03752 colvar%ss=s1(1)/s1(2)/REAL(nconf-1,dp) 03753 03754 ALLOCATE(fi(3,colvar%n_atom_s),stat=stat) 03755 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03756 03757 ii=0 03758 DO i=1,ncolv 03759 DO j=1,colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s 03760 ii=ii+1 03761 fi(:,ii)=colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:,j)*lambda*(ds1(i,1) & 03762 /s1(2)/REAL(nconf-1,dp)-colvar%ss*ds1(i,2)/s1(2))*2.0_dp 03763 END DO 03764 END DO 03765 03766 DO i=1,colvar%n_atom_s 03767 CALL put_derivative(colvar,i,fi(:,i)) 03768 END DO 03769 03770 DEALLOCATE(fi,stat=stat) 03771 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03772 DEALLOCATE(f_vals,stat=stat) 03773 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03774 DEALLOCATE(ss_vals,stat=stat) 03775 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03776 DEALLOCATE(s1v,stat=stat) 03777 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03778 DEALLOCATE(ds1v,stat=stat) 03779 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03780 DEALLOCATE(s1,stat=stat) 03781 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03782 DEALLOCATE(ds1,stat=stat) 03783 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03784 03785 END SUBROUTINE rpath_colvar 03786 03787 ! ***************************************************************************** 03794 SUBROUTINE rpath_dist_rmsd(colvar,particles,error) 03795 TYPE(colvar_type), POINTER :: colvar 03796 TYPE(particle_type), DIMENSION(:), 03797 POINTER :: particles 03798 TYPE(cp_error_type), INTENT(inout) :: error 03799 03800 CHARACTER(len=*), PARAMETER :: routineN = 'rpath_dist_rmsd', 03801 routineP = moduleN//':'//routineN 03802 03803 INTEGER :: i, iat, ii, ik, natom, nconf, 03804 rmsd_atom, stat 03805 INTEGER, DIMENSION(:), POINTER :: iatom 03806 LOGICAL :: failure 03807 REAL(dp) :: lambda, my_rmsd, s1(2), 03808 sum_exp 03809 REAL(dp), ALLOCATABLE, DIMENSION(:) :: r, r0, vec_dif 03810 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dvec_dif, fi, riat, s1v 03811 REAL(dp), ALLOCATABLE, 03812 DIMENSION(:, :, :) :: ds1 03813 REAL(dp), ALLOCATABLE, 03814 DIMENSION(:, :, :, :) :: ds1v 03815 REAL(dp), DIMENSION(:, :), POINTER :: path_conf 03816 03817 failure=.FALSE. 03818 03819 nconf=colvar%reaction_path_param%nr_frames 03820 rmsd_atom=colvar%reaction_path_param%n_components 03821 lambda=colvar%reaction_path_param%lambda 03822 path_conf => colvar%reaction_path_param%r_ref 03823 iatom => colvar%reaction_path_param%i_rmsd 03824 03825 natom = SIZE(particles) 03826 03827 ALLOCATE(r0(3*natom),STAT=stat) 03828 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03829 ALLOCATE(r(3*natom),STAT=stat) 03830 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03831 ALLOCATE(riat(3,rmsd_atom),STAT=stat) 03832 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03833 ALLOCATE(vec_dif(rmsd_atom),STAT=stat) 03834 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03835 ALLOCATE(dvec_dif(3,rmsd_atom),STAT=stat) 03836 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03837 ALLOCATE(s1v(2,nconf),STAT=stat) 03838 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03839 ALLOCATE(ds1v(3,rmsd_atom,2,nconf),STAT=stat) 03840 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03841 ALLOCATE(ds1(3,rmsd_atom,2),STAT=stat) 03842 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03843 DO i = 1,natom 03844 ii = (i-1)*3 03845 r0(ii+1) = particles(i)%r(1) 03846 r0(ii+2) = particles(i)%r(2) 03847 r0(ii+3) = particles(i)%r(3) 03848 END DO 03849 03850 DO iat = 1,rmsd_atom 03851 ii = iatom(iat) 03852 riat(:,iat) = particles(ii)%r 03853 END DO 03854 03855 DO ik = 1, nconf 03856 DO i = 1,natom 03857 ii = (i-1)*3 03858 r(ii+1) = path_conf(ii+1,ik) 03859 r(ii+2) = path_conf(ii+2,ik) 03860 r(ii+3) = path_conf(ii+3,ik) 03861 END DO 03862 03863 CALL rmsd3(particles,r,r0,output_unit=-1,my_val=my_rmsd,rotate=.TRUE.,error=error) 03864 03865 sum_exp = 0.0_dp 03866 DO iat = 1,rmsd_atom 03867 i = iatom(iat) 03868 ii = (i-1)*3 03869 vec_dif(iat) = (riat(1,iat)-r(ii+1))**2+(riat(2,iat)-r(ii+2))**2& 03870 +(riat(3,iat)-r(ii+3))**2 03871 sum_exp = sum_exp + vec_dif(iat) 03872 END DO 03873 03874 s1v(1,ik) = REAL(ik-1,dp)*EXP(-lambda*sum_exp) 03875 s1v(2,ik) = EXP(-lambda*sum_exp) 03876 DO iat = 1,rmsd_atom 03877 i = iatom(iat) 03878 ii = (i-1)*3 03879 ds1v(1,iat,1,ik) = r(ii+1)*s1v(1,ik) 03880 ds1v(1,iat,2,ik) = r(ii+1)*s1v(2,ik) 03881 ds1v(2,iat,1,ik) = r(ii+2)*s1v(1,ik) 03882 ds1v(2,iat,2,ik) = r(ii+2)*s1v(2,ik) 03883 ds1v(3,iat,1,ik) = r(ii+3)*s1v(1,ik) 03884 ds1v(3,iat,2,ik) = r(ii+3)*s1v(2,ik) 03885 END DO 03886 03887 END DO 03888 s1(1) = accurate_sum(s1v(1,:)) 03889 s1(2) = accurate_sum(s1v(2,:)) 03890 DO i = 1,2 03891 DO iat = 1,rmsd_atom 03892 ds1(1,iat,i) = accurate_sum(ds1v(1,iat,i,:)) 03893 ds1(2,iat,i) = accurate_sum(ds1v(2,iat,i,:)) 03894 ds1(3,iat,i) = accurate_sum(ds1v(3,iat,i,:)) 03895 END DO 03896 END DO 03897 03898 colvar%ss = s1(1)/s1(2)/REAL(nconf-1,dp) 03899 03900 ALLOCATE(fi(3,rmsd_atom),stat=stat) 03901 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03902 03903 DO iat = 1,rmsd_atom 03904 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)) 03905 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)) 03906 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)) 03907 CALL put_derivative(colvar,iat,fi(:,iat)) 03908 END DO 03909 03910 DEALLOCATE(fi,STAT=stat) 03911 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03912 DEALLOCATE(r0,STAT=stat) 03913 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03914 DEALLOCATE(r,STAT=stat) 03915 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03916 DEALLOCATE(riat,STAT=stat) 03917 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03918 DEALLOCATE(vec_dif,STAT=stat) 03919 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03920 DEALLOCATE(dvec_dif,STAT=stat) 03921 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03922 DEALLOCATE(s1v,STAT=stat) 03923 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03924 DEALLOCATE(ds1v,STAT=stat) 03925 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03926 DEALLOCATE(ds1,STAT=stat) 03927 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03928 03929 END SUBROUTINE rpath_dist_rmsd 03930 03931 03932 SUBROUTINE rpath_rmsd(colvar,particles,error) 03933 TYPE(colvar_type), POINTER :: colvar 03934 TYPE(particle_type), DIMENSION(:), 03935 POINTER :: particles 03936 TYPE(cp_error_type), INTENT(inout) :: error 03937 03938 CHARACTER(len=*), PARAMETER :: routineN = 'rpath_rmsd', 03939 routineP = moduleN//':'//routineN 03940 03941 INTEGER :: i, iat, ii, ik, natom, nconf, 03942 rmsd_atom, stat 03943 INTEGER, DIMENSION(:), POINTER :: iatom 03944 LOGICAL :: failure 03945 REAL(dp) :: lambda, my_rmsd, s1(2) 03946 REAL(dp), ALLOCATABLE, DIMENSION(:) :: r, r0 03947 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fi, riat, s1v 03948 REAL(dp), ALLOCATABLE, 03949 DIMENSION(:, :, :) :: ds1 03950 REAL(dp), ALLOCATABLE, 03951 DIMENSION(:, :, :, :) :: ds1v 03952 REAL(dp), DIMENSION(:, :), POINTER :: path_conf 03953 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weight 03954 REAL(KIND=dp), ALLOCATABLE, 03955 DIMENSION(:, :) :: drmsd 03956 03957 failure=.FALSE. 03958 03959 nconf=colvar%reaction_path_param%nr_frames 03960 rmsd_atom=colvar%reaction_path_param%n_components 03961 lambda=colvar%reaction_path_param%lambda 03962 path_conf => colvar%reaction_path_param%r_ref 03963 iatom => colvar%reaction_path_param%i_rmsd 03964 03965 natom = SIZE(particles) 03966 03967 ALLOCATE(r0(3*natom),STAT=stat) 03968 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03969 ALLOCATE(r(3*natom),STAT=stat) 03970 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03971 ALLOCATE(riat(3,rmsd_atom),STAT=stat) 03972 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03973 ALLOCATE(s1v(2,nconf),STAT=stat) 03974 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03975 ALLOCATE(ds1v(3,rmsd_atom,2,nconf),STAT=stat) 03976 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03977 ALLOCATE(ds1(3,rmsd_atom,2),STAT=stat) 03978 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03979 ALLOCATE(drmsd(3,natom),STAT=stat) 03980 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03981 drmsd = 0.0_dp 03982 ALLOCATE (weight(natom),STAT=stat) 03983 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 03984 03985 DO i = 1,natom 03986 ii = (i-1)*3 03987 r0(ii+1) = particles(i)%r(1) 03988 r0(ii+2) = particles(i)%r(2) 03989 r0(ii+3) = particles(i)%r(3) 03990 END DO 03991 03992 DO iat = 1,rmsd_atom 03993 ii = iatom(iat) 03994 riat(:,iat) = particles(ii)%r 03995 END DO 03996 03997 ! set weights of atoms in the rmsd list 03998 weight = 0.0_dp 03999 DO iat = 1,rmsd_atom 04000 i = iatom(iat) 04001 weight(i) = 1.0_dp 04002 END DO 04003 04004 04005 DO ik = 1, nconf 04006 DO i = 1,natom 04007 ii = (i-1)*3 04008 r(ii+1) = path_conf(ii+1,ik) 04009 r(ii+2) = path_conf(ii+2,ik) 04010 r(ii+3) = path_conf(ii+3,ik) 04011 END DO 04012 04013 CALL rmsd3(particles,r0,r,output_unit=-1,weights=weight,my_val=my_rmsd,& 04014 rotate=.FALSE.,drmsd3=drmsd,error=error) 04015 04016 s1v(1,ik) = REAL(ik-1,dp)*EXP(-lambda*my_rmsd) 04017 s1v(2,ik) = EXP(-lambda*my_rmsd) 04018 DO iat = 1,rmsd_atom 04019 i = iatom(iat) 04020 ds1v(1,iat,1,ik) = drmsd(1,i)*s1v(1,ik) 04021 ds1v(1,iat,2,ik) = drmsd(1,i)*s1v(2,ik) 04022 ds1v(2,iat,1,ik) = drmsd(2,i)*s1v(1,ik) 04023 ds1v(2,iat,2,ik) = drmsd(2,i)*s1v(2,ik) 04024 ds1v(3,iat,1,ik) = drmsd(3,i)*s1v(1,ik) 04025 ds1v(3,iat,2,ik) = drmsd(3,i)*s1v(2,ik) 04026 END DO 04027 END DO ! ik 04028 04029 s1(1) = accurate_sum(s1v(1,:)) 04030 s1(2) = accurate_sum(s1v(2,:)) 04031 DO i = 1,2 04032 DO iat = 1,rmsd_atom 04033 ds1(1,iat,i) = accurate_sum(ds1v(1,iat,i,:)) 04034 ds1(2,iat,i) = accurate_sum(ds1v(2,iat,i,:)) 04035 ds1(3,iat,i) = accurate_sum(ds1v(3,iat,i,:)) 04036 END DO 04037 END DO 04038 04039 colvar%ss = s1(1)/s1(2)/REAL(nconf-1,dp) 04040 04041 ALLOCATE(fi(3,rmsd_atom),stat=stat) 04042 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04043 04044 DO iat = 1,rmsd_atom 04045 fi(1,iat) = -lambda/s1(2)/REAL(nconf-1,dp)*(ds1(1,iat,1) - ds1(1,iat,2)*s1(1)/s1(2)) 04046 fi(2,iat) = -lambda/s1(2)/REAL(nconf-1,dp)*(ds1(2,iat,1) - ds1(2,iat,2)*s1(1)/s1(2)) 04047 fi(3,iat) = -lambda/s1(2)/REAL(nconf-1,dp)*(ds1(3,iat,1) - ds1(3,iat,2)*s1(1)/s1(2)) 04048 CALL put_derivative(colvar,iat,fi(:,iat)) 04049 END DO 04050 04051 DEALLOCATE(fi,STAT=stat) 04052 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04053 DEALLOCATE(r0,STAT=stat) 04054 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04055 DEALLOCATE(r,STAT=stat) 04056 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04057 DEALLOCATE(riat,STAT=stat) 04058 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04059 DEALLOCATE(s1v,STAT=stat) 04060 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04061 DEALLOCATE(ds1v,STAT=stat) 04062 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04063 DEALLOCATE(ds1,STAT=stat) 04064 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04065 DEALLOCATE(drmsd,STAT=stat) 04066 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04067 DEALLOCATE(weight,STAT=stat) 04068 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04069 04070 END SUBROUTINE rpath_rmsd 04071 04072 ! ***************************************************************************** 04078 SUBROUTINE distance_from_path_colvar(colvar,cell,subsys,particles,error) 04079 TYPE(colvar_type), POINTER :: colvar 04080 TYPE(cell_type), POINTER :: cell 04081 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04082 TYPE(particle_type), DIMENSION(:), 04083 OPTIONAL, POINTER :: particles 04084 TYPE(cp_error_type), INTENT(inout) :: error 04085 04086 CHARACTER(len=*), PARAMETER :: routineN = 'distance_from_path_colvar', 04087 routineP = moduleN//':'//routineN 04088 04089 LOGICAL :: failure 04090 TYPE(particle_list_type), POINTER :: particles_i 04091 TYPE(particle_type), DIMENSION(:), 04092 POINTER :: my_particles 04093 04094 failure=.FALSE. 04095 04096 CPPrecondition(colvar%type_id==distance_from_path_colvar_id,cp_failure_level,routineP,error,failure) 04097 IF (PRESENT(particles)) THEN 04098 my_particles => particles 04099 ELSE 04100 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 04101 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 04102 my_particles => particles_i%els 04103 END IF 04104 04105 IF(colvar%reaction_path_param%dist_rmsd) THEN 04106 CALL dpath_dist_rmsd(colvar,my_particles,error=error) 04107 ELSEIF(colvar%reaction_path_param%rmsd) THEN 04108 CALL dpath_rmsd(colvar,my_particles,error=error) 04109 ELSE 04110 CALL dpath_colvar(colvar,cell,my_particles,error=error) 04111 END IF 04112 04113 END SUBROUTINE distance_from_path_colvar 04114 04115 ! ***************************************************************************** 04122 SUBROUTINE dpath_colvar(colvar,cell,particles,error) 04123 TYPE(colvar_type), POINTER :: colvar 04124 TYPE(cell_type), POINTER :: cell 04125 TYPE(particle_type), DIMENSION(:), 04126 POINTER :: particles 04127 TYPE(cp_error_type), INTENT(inout) :: error 04128 04129 CHARACTER(len=*), PARAMETER :: routineN = 'dpath_colvar', 04130 routineP = moduleN//':'//routineN 04131 04132 INTEGER :: i, iend, ii, istart, j, k, 04133 ncolv, stat 04134 LOGICAL :: failure 04135 REAL(dp) :: lambda, s1 04136 REAL(dp), ALLOCATABLE, DIMENSION(:) :: ds1, s1v, ss_vals 04137 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ds1v, f_vals, fi 04138 04139 istart=colvar%reaction_path_param%function_bounds(1) 04140 iend=colvar%reaction_path_param%function_bounds(2) 04141 04142 ncolv=colvar%reaction_path_param%n_components 04143 lambda=colvar%reaction_path_param%lambda 04144 ALLOCATE(f_vals(ncolv,istart:iend),stat=stat) 04145 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04146 f_vals=colvar%reaction_path_param%f_vals 04147 ALLOCATE(ss_vals(ncolv),stat=stat) 04148 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04149 04150 DO i=1,ncolv 04151 CALL colvar_recursive_eval(colvar%reaction_path_param%colvar_p(i)%colvar,cell,particles,error) 04152 ss_vals(i)=colvar%reaction_path_param%colvar_p(i)%colvar%ss 04153 ENDDO 04154 04155 ALLOCATE(s1v(istart:iend),stat=stat) 04156 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04157 ALLOCATE(ds1v(ncolv,istart:iend),stat=stat) 04158 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04159 ALLOCATE(ds1(ncolv),stat=stat) 04160 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04161 04162 DO k=istart,iend 04163 s1v(k)=EXP(-lambda*DOT_PRODUCT(ss_vals(:)-f_vals(:,k),ss_vals(:)-f_vals(:,k))) 04164 DO j=1,ncolv 04165 ds1v(j,k)=f_vals(j,k)*s1v(k) 04166 END DO 04167 END DO 04168 04169 s1 = accurate_sum(s1v(:)) 04170 DO j=1,ncolv 04171 ds1(j) = accurate_sum(ds1v(j,:)) 04172 END DO 04173 colvar%ss=-1.0_dp/lambda*LOG(s1) 04174 04175 ALLOCATE(fi(3,colvar%n_atom_s),stat=stat) 04176 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04177 04178 ii=0 04179 DO i=1,ncolv 04180 DO j=1,colvar%reaction_path_param%colvar_p(i)%colvar%n_atom_s 04181 ii=ii+1 04182 fi(:,ii)=colvar%reaction_path_param%colvar_p(i)%colvar%dsdr(:,j)*& 04183 2.0_dp*( ss_vals(i) - ds1(i)/s1) 04184 END DO 04185 END DO 04186 04187 DO i=1,colvar%n_atom_s 04188 CALL put_derivative(colvar,i,fi(:,i)) 04189 END DO 04190 04191 DEALLOCATE(fi,stat=stat) 04192 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04193 DEALLOCATE(f_vals,stat=stat) 04194 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04195 DEALLOCATE(ss_vals,stat=stat) 04196 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04197 DEALLOCATE(s1v,stat=stat) 04198 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04199 DEALLOCATE(ds1v,stat=stat) 04200 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04201 DEALLOCATE(ds1,stat=stat) 04202 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04203 04204 END SUBROUTINE dpath_colvar 04205 04206 ! ***************************************************************************** 04213 SUBROUTINE dpath_dist_rmsd(colvar,particles,error) 04214 04215 TYPE(colvar_type), POINTER :: colvar 04216 TYPE(particle_type), DIMENSION(:), 04217 POINTER :: particles 04218 TYPE(cp_error_type), INTENT(inout) :: error 04219 04220 CHARACTER(len=*), PARAMETER :: routineN = 'dpath_dist_rmsd', 04221 routineP = moduleN//':'//routineN 04222 04223 INTEGER :: i, iat, ii, ik, natom, nconf, 04224 rmsd_atom, stat 04225 INTEGER, DIMENSION(:), POINTER :: iatom 04226 LOGICAL :: failure 04227 REAL(dp) :: lambda, s1, sum_exp 04228 REAL(dp), ALLOCATABLE, DIMENSION(:) :: r, r0, s1v, vec_dif 04229 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ds1, dvec_dif, fi, riat 04230 REAL(dp), ALLOCATABLE, 04231 DIMENSION(:, :, :) :: ds1v 04232 REAL(dp), DIMENSION(:, :), POINTER :: path_conf 04233 04234 failure=.FALSE. 04235 04236 nconf=colvar%reaction_path_param%nr_frames 04237 rmsd_atom=colvar%reaction_path_param%n_components 04238 lambda=colvar%reaction_path_param%lambda 04239 path_conf => colvar%reaction_path_param%r_ref 04240 iatom => colvar%reaction_path_param%i_rmsd 04241 04242 natom = SIZE(particles) 04243 04244 ALLOCATE(r0(3*natom),STAT=stat) 04245 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04246 ALLOCATE(r(3*natom),STAT=stat) 04247 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04248 ALLOCATE(riat(3,rmsd_atom),STAT=stat) 04249 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04250 ALLOCATE(vec_dif(rmsd_atom),STAT=stat) 04251 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04252 ALLOCATE(dvec_dif(3,rmsd_atom),STAT=stat) 04253 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04254 ALLOCATE(s1v(nconf),STAT=stat) 04255 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04256 ALLOCATE(ds1v(3,rmsd_atom,nconf),STAT=stat) 04257 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04258 ALLOCATE(ds1(3,rmsd_atom),STAT=stat) 04259 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04260 DO i = 1,natom 04261 ii = (i-1)*3 04262 r0(ii+1) = particles(i)%r(1) 04263 r0(ii+2) = particles(i)%r(2) 04264 r0(ii+3) = particles(i)%r(3) 04265 END DO 04266 04267 DO iat = 1,rmsd_atom 04268 ii = iatom(iat) 04269 riat(:,iat) = particles(ii)%r 04270 END DO 04271 04272 DO ik = 1, nconf 04273 DO i = 1,natom 04274 ii = (i-1)*3 04275 r(ii+1) = path_conf(ii+1,ik) 04276 r(ii+2) = path_conf(ii+2,ik) 04277 r(ii+3) = path_conf(ii+3,ik) 04278 END DO 04279 04280 CALL rmsd3(particles,r,r0,output_unit=-1,rotate=.TRUE.,error=error) 04281 04282 sum_exp = 0.0_dp 04283 DO iat = 1,rmsd_atom 04284 i = iatom(iat) 04285 ii = (i-1)*3 04286 vec_dif(iat) = (riat(1,iat)-r(ii+1))**2+(riat(2,iat)-r(ii+2))**2+(riat(3,iat)-r(ii+3))**2 04287 sum_exp = sum_exp + vec_dif(iat) 04288 dvec_dif(1,iat) = r(ii+1) 04289 dvec_dif(2,iat) = r(ii+2) 04290 dvec_dif(3,iat) = r(ii+3) 04291 END DO 04292 s1v(ik)=EXP(-lambda*sum_exp) 04293 DO iat = 1,rmsd_atom 04294 ds1v(1,iat,ik) = dvec_dif(1,iat) * s1v(ik) 04295 ds1v(2,iat,ik) = dvec_dif(2,iat) * s1v(ik) 04296 ds1v(3,iat,ik) = dvec_dif(3,iat) * s1v(ik) 04297 END DO 04298 END DO 04299 04300 s1 = accurate_sum(s1v(:)) 04301 DO iat = 1,rmsd_atom 04302 ds1(1,iat) = accurate_sum( ds1v(1,iat,:)) 04303 ds1(2,iat) = accurate_sum( ds1v(2,iat,:)) 04304 ds1(3,iat) = accurate_sum( ds1v(3,iat,:)) 04305 END DO 04306 colvar%ss=-1.0_dp/lambda*LOG(s1) 04307 04308 ALLOCATE(fi(3,rmsd_atom),stat=stat) 04309 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04310 04311 DO iat = 1,rmsd_atom 04312 fi(:,iat) = 2.0_dp*(riat(:,iat)-ds1(:,iat)/s1) 04313 CALL put_derivative(colvar,iat,fi(:,iat)) 04314 END DO 04315 04316 DEALLOCATE(fi,STAT=stat) 04317 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04318 DEALLOCATE(r0,STAT=stat) 04319 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04320 DEALLOCATE(r,STAT=stat) 04321 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04322 DEALLOCATE(riat,STAT=stat) 04323 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04324 DEALLOCATE(vec_dif,STAT=stat) 04325 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04326 DEALLOCATE(dvec_dif,STAT=stat) 04327 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04328 DEALLOCATE(s1v,STAT=stat) 04329 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04330 DEALLOCATE(ds1v,STAT=stat) 04331 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04332 DEALLOCATE(ds1,STAT=stat) 04333 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04334 END SUBROUTINE dpath_dist_rmsd 04335 04336 SUBROUTINE dpath_rmsd(colvar,particles,error) 04337 04338 TYPE(colvar_type), POINTER :: colvar 04339 TYPE(particle_type), DIMENSION(:), 04340 POINTER :: particles 04341 TYPE(cp_error_type), INTENT(inout) :: error 04342 04343 CHARACTER(len=*), PARAMETER :: routineN = 'dpath_rmsd', 04344 routineP = moduleN//':'//routineN 04345 04346 INTEGER :: i, iat, ii, ik, natom, nconf, 04347 rmsd_atom, stat 04348 INTEGER, DIMENSION(:), POINTER :: iatom 04349 LOGICAL :: failure 04350 REAL(dp) :: lambda, my_rmsd, s1 04351 REAL(dp), ALLOCATABLE, DIMENSION(:) :: r, r0, s1v 04352 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ds1, fi, riat 04353 REAL(dp), ALLOCATABLE, 04354 DIMENSION(:, :, :) :: ds1v 04355 REAL(dp), DIMENSION(:, :), POINTER :: path_conf 04356 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weight 04357 REAL(KIND=dp), ALLOCATABLE, 04358 DIMENSION(:, :) :: drmsd 04359 04360 failure=.FALSE. 04361 04362 nconf=colvar%reaction_path_param%nr_frames 04363 rmsd_atom=colvar%reaction_path_param%n_components 04364 lambda=colvar%reaction_path_param%lambda 04365 path_conf => colvar%reaction_path_param%r_ref 04366 iatom => colvar%reaction_path_param%i_rmsd 04367 04368 natom = SIZE(particles) 04369 04370 ALLOCATE(r0(3*natom),STAT=stat) 04371 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04372 ALLOCATE(r(3*natom),STAT=stat) 04373 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04374 ALLOCATE(riat(3,rmsd_atom),STAT=stat) 04375 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04376 ALLOCATE(s1v(nconf),STAT=stat) 04377 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04378 ALLOCATE(ds1v(3,rmsd_atom,nconf),STAT=stat) 04379 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04380 ALLOCATE(ds1(3,rmsd_atom),STAT=stat) 04381 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04382 ALLOCATE(drmsd(3,natom),STAT=stat) 04383 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04384 drmsd = 0.0_dp 04385 ALLOCATE (weight(natom),STAT=stat) 04386 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04387 04388 04389 DO i = 1,natom 04390 ii = (i-1)*3 04391 r0(ii+1) = particles(i)%r(1) 04392 r0(ii+2) = particles(i)%r(2) 04393 r0(ii+3) = particles(i)%r(3) 04394 END DO 04395 04396 DO iat = 1,rmsd_atom 04397 ii = iatom(iat) 04398 riat(:,iat) = particles(ii)%r 04399 END DO 04400 04401 ! set weights of atoms in the rmsd list 04402 weight = 0.0_dp 04403 DO iat = 1,rmsd_atom 04404 i = iatom(iat) 04405 weight(i) = 1.0_dp 04406 END DO 04407 04408 DO ik = 1, nconf 04409 DO i = 1,natom 04410 ii = (i-1)*3 04411 r(ii+1) = path_conf(ii+1,ik) 04412 r(ii+2) = path_conf(ii+2,ik) 04413 r(ii+3) = path_conf(ii+3,ik) 04414 END DO 04415 04416 CALL rmsd3(particles,r0,r,output_unit=-1,weights=weight,my_val=my_rmsd,& 04417 rotate=.FALSE.,drmsd3=drmsd,error=error) 04418 04419 s1v(ik)=EXP(-lambda*my_rmsd) 04420 DO iat = 1,rmsd_atom 04421 i = iatom(iat) 04422 ds1v(1,iat,ik) = drmsd(1,i) * s1v(ik) 04423 ds1v(2,iat,ik) = drmsd(2,i) * s1v(ik) 04424 ds1v(3,iat,ik) = drmsd(3,i) * s1v(ik) 04425 END DO 04426 END DO 04427 04428 s1 = accurate_sum(s1v(:)) 04429 DO iat = 1,rmsd_atom 04430 ds1(1,iat) = accurate_sum( ds1v(1,iat,:)) 04431 ds1(2,iat) = accurate_sum( ds1v(2,iat,:)) 04432 ds1(3,iat) = accurate_sum( ds1v(3,iat,:)) 04433 END DO 04434 colvar%ss=-1.0_dp/lambda*LOG(s1) 04435 04436 ALLOCATE(fi(3,rmsd_atom),stat=stat) 04437 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04438 04439 DO iat = 1,rmsd_atom 04440 fi(:,iat) = ds1(:,iat)/s1 04441 CALL put_derivative(colvar,iat,fi(:,iat)) 04442 END DO 04443 04444 DEALLOCATE(fi,STAT=stat) 04445 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04446 DEALLOCATE(r0,STAT=stat) 04447 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04448 DEALLOCATE(r,STAT=stat) 04449 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04450 DEALLOCATE(riat,STAT=stat) 04451 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04452 DEALLOCATE(s1v,STAT=stat) 04453 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04454 DEALLOCATE(ds1v,STAT=stat) 04455 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04456 DEALLOCATE(ds1,STAT=stat) 04457 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04458 DEALLOCATE(drmsd,STAT=stat) 04459 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04460 DEALLOCATE(weight,STAT=stat) 04461 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 04462 04463 END SUBROUTINE dpath_rmsd 04464 04465 ! ***************************************************************************** 04470 SUBROUTINE population_colvar(colvar,cell,subsys,particles,error) 04471 TYPE(colvar_type), POINTER :: colvar 04472 TYPE(cell_type), POINTER :: cell 04473 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04474 TYPE(particle_type), DIMENSION(:), 04475 OPTIONAL, POINTER :: particles 04476 TYPE(cp_error_type), INTENT(inout) :: error 04477 04478 CHARACTER(len=*), PARAMETER :: routineN = 'population_colvar', 04479 routineP = moduleN//':'//routineN 04480 04481 INTEGER :: i, ii, jj, n_atoms_from, 04482 n_atoms_to, ndcrd, nncrd 04483 LOGICAL :: failure 04484 REAL(dp) :: dfunc, dfunc_coord, ftmp(3), func, func_coord, 04485 inv_n_atoms_from, invden, n_0, ncoord, norm, num, population, r12, r_0, 04486 rdist, sigma, ss(3), xij(3) 04487 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftmp_coord 04488 REAL(dp), DIMENSION(3) :: xpi, xpj 04489 TYPE(particle_list_type), POINTER :: particles_i 04490 TYPE(particle_type), DIMENSION(:), 04491 POINTER :: my_particles 04492 04493 failure=.FALSE. 04494 ! If we defined the coordination number with KINDS then we have still 04495 ! to fill few missing informations... 04496 NULLIFY(particles_i) 04497 CPPrecondition(colvar%type_id==population_colvar_id,cp_failure_level,routineP,error,failure) 04498 IF (PRESENT(particles)) THEN 04499 my_particles => particles 04500 ELSE 04501 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 04502 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 04503 my_particles => particles_i%els 04504 END IF 04505 n_atoms_to=colvar%population_param%n_atoms_to 04506 n_atoms_from=colvar%population_param%n_atoms_from 04507 nncrd=colvar%population_param%nncrd 04508 ndcrd=colvar%population_param%ndcrd 04509 r_0=colvar%population_param%r_0 04510 n_0=colvar%population_param%n0 04511 sigma=colvar%population_param%sigma 04512 04513 ALLOCATE(ftmp_coord(3,n_atoms_to)) 04514 ftmp_coord=0.0_dp 04515 04516 ncoord=0.0_dp 04517 population=0.0_dp 04518 04519 colvar%dsdr=0.0_dp 04520 inv_n_atoms_from=1.0_dp/REAL(n_atoms_from) 04521 04522 norm=SQRT(pi*2.0_dp)*sigma 04523 norm=1/norm 04524 04525 DO ii=1,n_atoms_from 04526 i=colvar%population_param%i_at_from(ii) 04527 CALL get_coordinates(colvar, i, xpi, my_particles) 04528 DO jj=1,n_atoms_to 04529 i=colvar%population_param%i_at_to(jj) 04530 CALL get_coordinates(colvar, i, xpj, my_particles) 04531 ss=MATMUL(cell%h_inv,xpi(:)-xpj(:)) 04532 ss=ss-NINT(ss) 04533 xij=MATMUL(cell%hmat,ss) 04534 r12=SQRT(xij(1)**2+xij(2)**2+xij(3)**2) 04535 IF(r12 < 1.0e-8_dp)CYCLE 04536 rdist = r12/r_0 04537 num=(1.0_dp-rdist**nncrd) 04538 invden=1.0_dp/(1.0_dp-rdist**ndcrd) 04539 func_coord=num*invden 04540 dfunc_coord= (- nncrd*rdist **(nncrd-1) *invden & 04541 + num*(invden)**2 * ndcrd *rdist **(ndcrd-1))/(r12*r_0) 04542 04543 ncoord=ncoord+func_coord 04544 ftmp_coord(1,jj) = dfunc_coord*xij(1) 04545 ftmp_coord(2,jj) = dfunc_coord*xij(2) 04546 ftmp_coord(3,jj) = dfunc_coord*xij(3) 04547 END DO 04548 04549 func=EXP(-(ncoord-n_0)**2/(2.0_dp*sigma*sigma)) 04550 dfunc=-func*(ncoord-n_0)/(sigma*sigma) 04551 04552 population=population+norm*func 04553 DO jj=1,n_atoms_to 04554 ftmp(1)=ftmp_coord(1,jj)*dfunc 04555 ftmp(2)=ftmp_coord(2,jj)*dfunc 04556 ftmp(3)=ftmp_coord(3,jj)*dfunc 04557 CALL put_derivative(colvar, ii, ftmp) 04558 ftmp(1)=-ftmp_coord(1,jj)*dfunc 04559 ftmp(2)=-ftmp_coord(2,jj)*dfunc 04560 ftmp(3)=-ftmp_coord(3,jj)*dfunc 04561 CALL put_derivative(colvar, n_atoms_from+jj, ftmp) 04562 ENDDO 04563 ncoord=0.0_dp 04564 ENDDO 04565 colvar%ss=population 04566 END SUBROUTINE population_colvar 04567 04568 ! ***************************************************************************** 04574 SUBROUTINE gyration_radius_colvar(colvar,cell,subsys,particles,error) 04575 04576 TYPE(colvar_type), POINTER :: colvar 04577 TYPE(cell_type), POINTER :: cell 04578 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04579 TYPE(particle_type), DIMENSION(:), 04580 OPTIONAL, POINTER :: particles 04581 TYPE(cp_error_type), INTENT(inout) :: error 04582 04583 CHARACTER(len=*), PARAMETER :: routineN = 'gyration_radius_colvar', 04584 routineP = moduleN//':'//routineN 04585 04586 INTEGER :: i, ii, n_atoms 04587 LOGICAL :: failure 04588 REAL(dp) :: dri2, func, gyration, inv_n, 04589 mass_tot, mi 04590 REAL(dp), DIMENSION(3) :: dfunc, dxi, ftmp, ss, xpcom, 04591 xpi 04592 TYPE(particle_list_type), POINTER :: particles_i 04593 TYPE(particle_type), DIMENSION(:), 04594 POINTER :: my_particles 04595 04596 failure=.FALSE. 04597 04598 NULLIFY(particles_i,my_particles) 04599 CPPrecondition(colvar%type_id==gyration_colvar_id,cp_failure_level,routineP,error,failure) 04600 IF (PRESENT(particles)) THEN 04601 my_particles => particles 04602 ELSE 04603 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 04604 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 04605 my_particles => particles_i%els 04606 END IF 04607 n_atoms=colvar%gyration_param%n_atoms 04608 inv_n = 1.0_dp/n_atoms 04609 04610 04611 !compute COM position 04612 xpcom = 0.0_dp 04613 mass_tot = 0.0_dp 04614 DO ii = 1,n_atoms 04615 i=colvar%gyration_param%i_at(ii) 04616 CALL get_coordinates(colvar, i, xpi, my_particles) 04617 CALL get_mass(colvar, i, mi, my_particles) 04618 xpcom(:) = xpcom(:) + xpi(:)*mi 04619 mass_tot = mass_tot + mi 04620 END DO 04621 xpcom(:) = xpcom(:)/mass_tot 04622 04623 func = 0.0_dp 04624 ftmp = 0.0_dp 04625 dfunc = 0.0_dp 04626 DO ii = 1,n_atoms 04627 i=colvar%gyration_param%i_at(ii) 04628 CALL get_coordinates(colvar, i, xpi, my_particles) 04629 ss=MATMUL(cell%h_inv,xpi(:)-xpcom(:)) 04630 ss=ss-NINT(ss) 04631 dxi=MATMUL(cell%hmat,ss) 04632 dri2=(dxi(1)**2+dxi(2)**2+dxi(3)**2) 04633 func = func + dri2 04634 dfunc(:) = dfunc(:) + dxi(:) 04635 END DO 04636 gyration = SQRT(inv_n*func) 04637 04638 DO ii = 1,n_atoms 04639 i=colvar%gyration_param%i_at(ii) 04640 CALL get_coordinates(colvar, i, xpi, my_particles) 04641 CALL get_mass(colvar, i, mi, my_particles) 04642 ss=MATMUL(cell%h_inv,xpi(:)-xpcom(:)) 04643 ss=ss-NINT(ss) 04644 dxi=MATMUL(cell%hmat,ss) 04645 ftmp(1) = dxi(1) - dfunc(1)*mi/mass_tot 04646 ftmp(2) = dxi(2) - dfunc(2)*mi/mass_tot 04647 ftmp(3) = dxi(3) - dfunc(3)*mi/mass_tot 04648 ftmp(:) = ftmp(:)*inv_n/gyration 04649 CALL put_derivative(colvar, ii, ftmp) 04650 END DO 04651 colvar%ss=gyration 04652 04653 END SUBROUTINE gyration_radius_colvar 04654 04655 ! ***************************************************************************** 04661 SUBROUTINE rmsd_colvar(colvar,subsys,particles,error) 04662 TYPE(colvar_type), POINTER :: colvar 04663 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04664 TYPE(particle_type), DIMENSION(:), 04665 OPTIONAL, POINTER :: particles 04666 TYPE(cp_error_type), INTENT(inout) :: error 04667 04668 CHARACTER(len=*), PARAMETER :: routineN = 'rmsd_colvar', 04669 routineP = moduleN//':'//routineN 04670 04671 LOGICAL :: failure 04672 04673 failure = .FALSE. 04674 CALL rmsd_colvar_low(colvar,subsys,particles,error) 04675 END SUBROUTINE rmsd_colvar 04676 04677 ! ***************************************************************************** 04687 SUBROUTINE rmsd_colvar_low(colvar,subsys,particles,error) 04688 04689 TYPE(colvar_type), POINTER :: colvar 04690 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04691 TYPE(particle_type), DIMENSION(:), 04692 OPTIONAL, POINTER :: particles 04693 TYPE(cp_error_type), INTENT(inout) :: error 04694 04695 CHARACTER(len=*), PARAMETER :: routineN = 'rmsd_colvar_low', 04696 routineP = moduleN//':'//routineN 04697 04698 INTEGER :: i, ii, natom, nframes, stat 04699 LOGICAL :: failure 04700 REAL(kind=dp) :: cv_val, f1, ftmp(3) 04701 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: der, r, rmsd 04702 REAL(kind=dp), ALLOCATABLE, 04703 DIMENSION(:, :) :: r0 04704 REAL(kind=dp), ALLOCATABLE, 04705 DIMENSION(:, :, :) :: drmsd 04706 REAL(kind=dp), DIMENSION(:), POINTER :: weights 04707 TYPE(particle_list_type), POINTER :: particles_i 04708 TYPE(particle_type), DIMENSION(:), 04709 POINTER :: my_particles 04710 04711 failure=.FALSE. 04712 04713 NULLIFY(my_particles,particles_i,weights) 04714 CPPrecondition(colvar%type_id==rmsd_colvar_id,cp_failure_level,routineP,error,failure) 04715 IF (PRESENT(particles)) THEN 04716 my_particles => particles 04717 ELSE 04718 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 04719 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 04720 my_particles => particles_i%els 04721 END IF 04722 04723 natom = SIZE(my_particles) 04724 nframes=colvar%rmsd_param%nr_frames 04725 ALLOCATE (drmsd(3,natom,nframes),STAT=stat) 04726 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04727 drmsd = 0.0_dp 04728 04729 ALLOCATE (r0(3*natom,nframes),STAT=stat) 04730 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04731 ALLOCATE (rmsd(nframes),STAT=stat) 04732 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04733 ALLOCATE (der(nframes),STAT=stat) 04734 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04735 ALLOCATE (r(3*natom),STAT=stat) 04736 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04737 04738 weights => colvar%rmsd_param%weights 04739 DO i = 1,natom 04740 ii = (i-1)*3 04741 r(ii+1) = my_particles(i)%r(1) 04742 r(ii+2) = my_particles(i)%r(2) 04743 r(ii+3) = my_particles(i)%r(3) 04744 END DO 04745 r0 = colvar%rmsd_param%r_ref 04746 rmsd = 0.0_dp 04747 04748 CALL rmsd3(my_particles,r,r0(:,1),output_unit=-1,weights=weights,my_val=rmsd(1),rotate=.FALSE.,drmsd3=drmsd(:,:,1),error=error) 04749 04750 IF (nframes==2) THEN 04751 CALL rmsd3(my_particles,r,r0(:,2),output_unit=-1,weights=weights,my_val=rmsd(2),rotate=.FALSE.,drmsd3=drmsd(:,:,2),& 04752 error=error) 04753 04754 f1 = 1.0_dp/(rmsd(1)+rmsd(2)) 04755 ! (rmsdA-rmsdB)/(rmsdA+rmsdB) 04756 cv_val = (rmsd(1)-rmsd(2))*f1 04757 ! (rmsdA+rmsdB)^-1-(rmsdA-rmsdB)/(rmsdA+rmsdB)^2 04758 der(1) = f1-cv_val*f1 04759 ! -(rmsdA+rmsdB)^-1-(rmsdA-rmsdB)/(rmsdA+rmsdB)^2 04760 der(2) = -f1-cv_val*f1 04761 04762 DO i = 1,colvar%rmsd_param%n_atoms 04763 ii = colvar%rmsd_param%i_rmsd(i) 04764 IF(weights(ii)>0.0_dp) THEN 04765 ftmp(1) = der(1)*drmsd(1,ii,1)+der(2)*drmsd(1,ii,2) 04766 ftmp(2) = der(1)*drmsd(2,ii,1)+der(2)*drmsd(2,ii,2) 04767 ftmp(3) = der(1)*drmsd(3,ii,1)+der(2)*drmsd(3,ii,2) 04768 CALL put_derivative(colvar, i, ftmp) 04769 END IF 04770 END DO 04771 ELSE IF (nframes==1) THEN 04772 ! Protect in case of numerical issues (for two identical frames!) 04773 rmsd(1)=ABS(rmsd(1)) 04774 cv_val = SQRT(rmsd(1)) 04775 f1 = 0.0_dp 04776 IF (cv_val /= 0.0_dp) f1 = 0.5_dp/cv_val 04777 DO i = 1,colvar%rmsd_param%n_atoms 04778 ii = colvar%rmsd_param%i_rmsd(i) 04779 IF(weights(ii)>0.0_dp) THEN 04780 ftmp(1) = f1*drmsd(1,ii,1) 04781 ftmp(2) = f1*drmsd(2,ii,1) 04782 ftmp(3) = f1*drmsd(3,ii,1) 04783 CALL put_derivative(colvar, i, ftmp) 04784 END IF 04785 END DO 04786 ELSE 04787 CALL cp_unimplemented_error(fromWhere=routineP, & 04788 message="RMSD implemented only for 1 and 2 reference frames!", error=error, error_level=cp_failure_level) 04789 END IF 04790 colvar%ss=cv_val 04791 04792 DEALLOCATE(der,STAT=stat) 04793 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04794 DEALLOCATE(r0,STAT=stat) 04795 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04796 DEALLOCATE(r,STAT=stat) 04797 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04798 DEALLOCATE(drmsd,STAT=stat) 04799 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04800 DEALLOCATE(rmsd,STAT=stat) 04801 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 04802 04803 END SUBROUTINE rmsd_colvar_low 04804 04805 ! ***************************************************************************** 04811 SUBROUTINE ring_puckering_colvar(colvar,cell,subsys,particles,error) 04812 TYPE(colvar_type), POINTER :: colvar 04813 TYPE(cell_type), POINTER :: cell 04814 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 04815 TYPE(particle_type), DIMENSION(:), 04816 OPTIONAL, POINTER :: particles 04817 TYPE(cp_error_type), INTENT(inout) :: error 04818 04819 CHARACTER(len=*), PARAMETER :: routineN = 'ring_puckering_colvar', 04820 routineP = moduleN//':'//routineN 04821 04822 INTEGER :: i, ii, istat, j, jj, m, nring 04823 LOGICAL :: failure 04824 REAL(KIND=dp) :: a, at, b, da, db, ds, kr, 04825 rpxpp, svar 04826 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosj, sinj, z 04827 REAL(KIND=dp), ALLOCATABLE, 04828 DIMENSION(:, :) :: r 04829 REAL(KIND=dp), ALLOCATABLE, 04830 DIMENSION(:, :, :) :: nforce, zforce 04831 REAL(KIND=dp), DIMENSION(3) :: ftmp, nv, r0, rp, rpp, uv 04832 REAL(KIND=dp), DIMENSION(3, 3) :: dnvp, dnvpp 04833 TYPE(particle_list_type), POINTER :: particles_i 04834 TYPE(particle_type), DIMENSION(:), 04835 POINTER :: my_particles 04836 04837 failure=.FALSE. 04838 04839 CPPrecondition(colvar%type_id==ring_puckering_colvar_id,cp_failure_level,routineP,error,failure) 04840 IF (PRESENT(particles)) THEN 04841 my_particles => particles 04842 ELSE 04843 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 04844 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 04845 my_particles => particles_i%els 04846 END IF 04847 04848 nring = colvar%ring_puckering_param%nring 04849 ALLOCATE (r(3,nring),z(nring),cosj(nring),sinj(nring),STAT=istat) 04850 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 04851 ALLOCATE (nforce(3,3,nring),zforce(nring,nring,3),STAT=istat) 04852 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 04853 DO ii = 1,nring 04854 i=colvar%ring_puckering_param%atoms(ii) 04855 CALL get_coordinates(colvar, i, r(:,ii), my_particles) 04856 END DO 04857 ! get all atoms within PBC distance of atom 1 04858 r0(:) = r(:,1) 04859 DO ii = 1,nring 04860 r(:,ii) = pbc(r(:,ii),r0,cell) 04861 END DO 04862 !compute origin position 04863 r0 = 0.0_dp 04864 DO ii = 1,nring 04865 r0(:) = r0(:) + r(:,ii) 04866 END DO 04867 kr = 1._dp/REAL(nring,KIND=dp) 04868 r0(:) = r0(:)*kr 04869 DO ii = 1,nring 04870 r(:,ii) = r(:,ii) - r0(:) 04871 END DO 04872 ! orientation vectors 04873 rp = 0._dp 04874 rpp = 0._dp 04875 DO ii = 1,nring 04876 cosj(ii) = COS(twopi*(ii-1)*kr) 04877 sinj(ii) = SIN(twopi*(ii-1)*kr) 04878 rp(:) = rp(:) + r(:,ii)*sinj(ii) 04879 rpp(:) = rpp(:) + r(:,ii)*cosj(ii) 04880 END DO 04881 nv = vector_product(rp,rpp) 04882 nv = nv/SQRT(SUM(nv**2)) 04883 04884 ! derivatives of normal 04885 uv = vector_product(rp,rpp) 04886 rpxpp = SQRT(SUM(uv**2)) 04887 DO i=1,3 04888 uv = 0._dp 04889 uv(i) = 1._dp 04890 uv = vector_product(uv,rpp)/rpxpp 04891 dnvp(:,i) = uv - nv * SUM(uv*nv) 04892 uv = 0._dp 04893 uv(i) = 1._dp 04894 uv = vector_product(rp,uv)/rpxpp 04895 dnvpp(:,i) = uv - nv * SUM(uv*nv) 04896 END DO 04897 DO ii = 1,nring 04898 nforce(:,:,ii) = dnvp(:,:)*sinj(ii) + dnvpp(:,:)*cosj(ii) 04899 END DO 04900 04901 ! molecular z-coordinate 04902 DO ii = 1,nring 04903 z(ii) = SUM(r(:,ii)*nv(:)) 04904 END DO 04905 ! z-force 04906 DO ii = 1,nring 04907 DO jj = 1,nring 04908 IF(ii==jj) THEN 04909 zforce(ii,jj,:) = nv 04910 ELSE 04911 zforce(ii,jj,:) = 0._dp 04912 END IF 04913 DO i=1,3 04914 DO j=1,3 04915 zforce(ii,jj,i) = zforce(ii,jj,i) + r(j,ii)*nforce(j,i,jj) 04916 END DO 04917 END DO 04918 END DO 04919 END DO 04920 04921 IF(colvar%ring_puckering_param%iq == 0) THEN 04922 ! total puckering amplitude 04923 svar = SQRT(SUM(z**2)) 04924 DO ii = 1,nring 04925 ftmp = 0._dp 04926 DO jj = 1,nring 04927 ftmp(:) = ftmp(:) + zforce(jj,ii,:)*z(jj) 04928 END DO 04929 ftmp = ftmp/svar 04930 CALL put_derivative(colvar, ii, ftmp) 04931 END DO 04932 ELSE 04933 m = ABS(colvar%ring_puckering_param%iq) 04934 CPPrecondition(m/=1,cp_failure_level,routineP,error,failure) 04935 IF(MOD(nring,2)==0 .AND. colvar%ring_puckering_param%iq == nring/2) THEN 04936 ! single puckering amplitude 04937 svar = 0._dp 04938 DO ii = 1,nring 04939 IF(MOD(ii,2)==0) THEN 04940 svar = svar - z(ii) 04941 ELSE 04942 svar = svar + z(ii) 04943 END IF 04944 END DO 04945 svar = svar*SQRT(kr) 04946 DO ii = 1,nring 04947 ftmp = 0._dp 04948 DO jj=1,nring 04949 IF(MOD(jj,2)==0) THEN 04950 ftmp(:) = ftmp(:) - zforce(jj,ii,:)*SQRT(kr) 04951 ELSE 04952 ftmp(:) = ftmp(:) + zforce(jj,ii,:)*SQRT(kr) 04953 END IF 04954 END DO 04955 CALL put_derivative(colvar, ii, -ftmp) 04956 END DO 04957 ELSE 04958 CPPrecondition(m<=(nring-1)/2,cp_failure_level,routineP,error,failure) 04959 a = 0._dp 04960 b = 0._dp 04961 DO ii = 1,nring 04962 a = a + z(ii)*COS(twopi*m*(ii-1)*kr) 04963 b = b - z(ii)*SIN(twopi*m*(ii-1)*kr) 04964 END DO 04965 a = a*SQRT(2._dp*kr) 04966 b = b*SQRT(2._dp*kr) 04967 IF(colvar%ring_puckering_param%iq > 0) THEN 04968 ! puckering amplitude 04969 svar = SQRT(a*a+b*b) 04970 da = a/svar 04971 db = b/svar 04972 ELSE 04973 ! puckering phase angle 04974 at = ATAN2(a,b) 04975 IF(at>pi/2._dp) THEN 04976 svar = 2.5_dp*pi - at 04977 ELSE 04978 svar = 0.5_dp*pi - at 04979 END IF 04980 da = -b/(a*a+b*b) 04981 db = a/(a*a+b*b) 04982 END IF 04983 DO jj = 1,nring 04984 ftmp = 0._dp 04985 DO ii = 1,nring 04986 ds = da*COS(twopi*m*(ii-1)*kr) 04987 ds = ds - db*SIN(twopi*m*(ii-1)*kr) 04988 ftmp(:) = ftmp(:) + ds*SQRT(2._dp*kr)*zforce(ii,jj,:) 04989 END DO 04990 CALL put_derivative(colvar, jj, ftmp) 04991 END DO 04992 END IF 04993 END IF 04994 04995 colvar%ss=svar 04996 04997 DEALLOCATE (r,z,cosj,sinj,nforce,zforce,STAT=istat) 04998 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 04999 05000 END SUBROUTINE ring_puckering_colvar 05001 05002 ! ***************************************************************************** 05006 RECURSIVE FUNCTION rec_eval_grid(iw1,ncol,f_vals,v_count,& 05007 gp,grid_sp,step_size,istart,iend,s1v,s1,p_bounds,lambda,ifunc,nconf) RESULT(k) 05008 INTEGER :: iw1, ncol 05009 REAL(dp), DIMENSION(:, :), POINTER :: f_vals 05010 INTEGER :: v_count 05011 REAL(dp), DIMENSION(:), POINTER :: gp, grid_sp 05012 REAL(dp) :: step_size 05013 INTEGER :: istart, iend 05014 REAL(dp), DIMENSION(:, :), POINTER :: s1v 05015 REAL(dp), DIMENSION(:), POINTER :: s1 05016 INTEGER, DIMENSION(:, :), POINTER :: p_bounds 05017 REAL(dp) :: lambda 05018 INTEGER :: ifunc, nconf, k 05019 05020 INTEGER :: count1, i 05021 05022 k=1 05023 IF(v_count.LT.ncol)THEN 05024 count1=v_count+1 05025 DO i=p_bounds(1,count1),p_bounds(2,count1) 05026 gp(count1)=REAL(i,KIND=dp)*grid_sp(count1) 05027 k=rec_eval_grid(iw1,ncol,f_vals,count1,gp,grid_sp,step_size,istart,iend,s1v,s1,p_bounds,lambda,ifunc,nconf) 05028 END DO 05029 ELSE IF(v_count==ncol .AND. ifunc==1)THEN 05030 DO i=istart,iend 05031 s1v(1,i)=REAL(i,kind=dp)*step_size*EXP(-lambda*DOT_PRODUCT(gp(:)-f_vals(:,i), 05032 gp(:)-f_vals(:,i))) 05033 s1v(2,i)=EXP(-lambda*DOT_PRODUCT(gp(:)-f_vals(:,i),gp(:)-f_vals(:,i))) 05034 END DO 05035 DO i=1,2 05036 s1(i)=accurate_sum(s1v(i,:)) 05037 END DO 05038 WRITE(iw1, '(5F10.5)')gp(:),s1(1)/s1(2)/REAL(nconf-1,dp) 05039 ELSE IF(v_count==ncol .AND. ifunc==2)THEN 05040 DO i=istart,iend 05041 s1v(1,i)=EXP(-lambda*DOT_PRODUCT(gp(:)-f_vals(:,i),gp(:)-f_vals(:,i))) 05042 END DO 05043 s1(1)=accurate_sum(s1v(1,:)) 05044 05045 WRITE(iw1, '(5F10.5)')gp(:),-lambda*LOG(s1(1)) 05046 END IF 05047 END FUNCTION rec_eval_grid 05048 05049 ! ***************************************************************************** 05055 SUBROUTINE read_frames(frame_section,para_env,nr_frames,r_ref,n_atoms,error) 05056 05057 TYPE(section_vals_type), POINTER :: frame_section 05058 TYPE(cp_para_env_type), POINTER :: para_env 05059 INTEGER, INTENT(IN) :: nr_frames 05060 REAL(dp), DIMENSION(:, :), POINTER :: r_ref 05061 INTEGER, INTENT(OUT) :: n_atoms 05062 TYPE(cp_error_type), INTENT(inout) :: error 05063 05064 CHARACTER(len=*), PARAMETER :: routineN = 'read_frames', 05065 routineP = moduleN//':'//routineN 05066 05067 CHARACTER(LEN=default_path_length) :: filename 05068 CHARACTER(LEN=default_string_length) :: dummy_char 05069 INTEGER :: i, j, natom, stat 05070 LOGICAL :: explicit, failure, my_end 05071 REAL(KIND=dp), DIMENSION(:), POINTER :: rptr 05072 TYPE(cp_parser_type), POINTER :: parser 05073 TYPE(section_vals_type), POINTER :: coord_section 05074 05075 failure = .FALSE. 05076 NULLIFY(rptr) 05077 05078 DO i = 1,nr_frames 05079 coord_section => section_vals_get_subs_vals(frame_section,"COORD",i_rep_section=i,error=error) 05080 CALL section_vals_get(coord_section,explicit=explicit,error=error) 05081 ! Cartesian Coordinates 05082 IF (explicit) THEN 05083 CALL section_vals_val_get(coord_section,"_DEFAULT_KEYWORD_",& 05084 n_rep_val=natom,error=error) 05085 IF(i==1) THEN 05086 ALLOCATE(r_ref(3*natom,nr_frames), STAT=stat) 05087 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 05088 n_atoms = natom 05089 ELSE 05090 CPPostcondition(3*natom==SIZE(r_ref,1),cp_failure_level,routineP,error,failure) 05091 END IF 05092 DO j = 1, natom 05093 CALL section_vals_val_get(coord_section,"_DEFAULT_KEYWORD_",& 05094 i_rep_val=j,r_vals=rptr,error=error) 05095 r_ref((j-1)*3+1:(j-1)*3+3,i) = rptr(1:3) 05096 END DO ! natom 05097 ELSE 05098 CALL section_vals_val_get(frame_section,"COORD_FILE_NAME",i_rep_section=i,c_val=filename,error=error) 05099 CPPostcondition(TRIM(filename)/="",cp_failure_level,routineP,error,failure) 05100 NULLIFY(parser) 05101 ALLOCATE(rptr(3)) 05102 CALL parser_create(parser,filename,para_env=para_env,parse_white_lines=.TRUE.,& 05103 error=error) 05104 CALL parser_get_next_line(parser,1,error=error) 05105 ! Start parser 05106 CALL parser_get_object(parser,natom,error=error) 05107 CALL parser_get_next_line(parser,1,error=error) 05108 IF(i==1) THEN 05109 ALLOCATE(r_ref(3*natom,nr_frames), STAT=stat) 05110 CPPostcondition(stat == 0,cp_failure_level,routinep,error,failure) 05111 n_atoms = natom 05112 ELSE 05113 CPPostcondition(3*natom==SIZE(r_ref,1),cp_failure_level,routineP,error,failure) 05114 END IF 05115 DO j = 1, natom 05116 ! Atom coordinates 05117 CALL parser_get_next_line(parser,1,at_end=my_end,error=error) 05118 CALL cp_assert(.NOT.my_end,cp_fatal_level,cp_assertion_failed,routineP,& 05119 "Number of lines in XYZ format not equal to the number of atoms."//& 05120 " Error in XYZ format for COORD_A (CV rmsd). Very probably the"//& 05121 " line with title is missing or is empty. Please check the XYZ file and rerun your job!"//& 05122 CPSourceFileRef,& 05123 only_ionode=.TRUE.) 05124 READ(parser%input_line,*) dummy_char,rptr(1:3) 05125 r_ref((j-1)*3+1,i) = cp_unit_to_cp2k(rptr(1),"angstrom",error=error) 05126 r_ref((j-1)*3+2,i) = cp_unit_to_cp2k(rptr(2),"angstrom",error=error) 05127 r_ref((j-1)*3+3,i) = cp_unit_to_cp2k(rptr(3),"angstrom",error=error) 05128 END DO ! natom 05129 CALL parser_release(parser,error=error) 05130 DEALLOCATE(rptr) 05131 END IF 05132 END DO ! nr_frames 05133 05134 END SUBROUTINE read_frames 05135 05136 05137 ! ***************************************************************************** 05141 SUBROUTINE Wc_colvar(colvar,cell,subsys,particles,qs_env,error) 05142 TYPE(colvar_type), POINTER :: colvar 05143 TYPE(cell_type), POINTER :: cell 05144 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 05145 TYPE(particle_type), DIMENSION(:), 05146 OPTIONAL, POINTER :: particles 05147 TYPE(qs_environment_type), POINTER, OPTIONAL :: qs_env ! optional just because I am lazy... but I should get rid of it... 05148 TYPE(cp_error_type), INTENT(inout) :: error 05149 05150 CHARACTER(len=*), PARAMETER :: routineN = 'Wc_colvar', 05151 routineP = moduleN//':'//routineN 05152 05153 INTEGER :: Od, H, Oa 05154 LOGICAL :: failure 05155 REAL(dp) :: rOd(3), rOa(3), rH(3), 05156 x,y,s(3),xv(3),dmin,amin,aux 05157 INTEGER :: idmin, iamin,i,j,ai 05158 TYPE(particle_list_type), POINTER :: particles_i 05159 TYPE(particle_type), DIMENSION(:), 05160 POINTER :: my_particles 05161 TYPE(wannier_centres_type), DIMENSION(:),POINTER :: wc 05162 INTEGER, ALLOCATABLE :: wcai(:),wcdi(:) ! contains the indeces of the wannier centres closed to the donor and acceptor 05163 INTEGER :: nwca,nwcd,ierror 05164 REAL(dp) :: rcut 05165 05166 failure=.FALSE. 05167 NULLIFY(particles_i,wc) 05168 05169 CPPrecondition(colvar%type_id==Wc_colvar_id,cp_failure_level,routineP,error,failure) 05170 IF (PRESENT(particles)) THEN 05171 my_particles => particles 05172 ELSE 05173 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 05174 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 05175 my_particles => particles_i%els 05176 END IF 05177 CALL get_qs_env(qs_env,WannierCentres=wc,error=error) 05178 rcut =colvar%Wc%rcut ! distances are in bohr as far as I remember 05179 Od=colvar%Wc%ids(1) 05180 H=colvar%Wc%ids(2) 05181 Oa=colvar%Wc%ids(3) 05182 CALL get_coordinates(colvar, Od, rOd, my_particles) 05183 CALL get_coordinates(colvar, H, rH, my_particles) 05184 CALL get_coordinates(colvar, Oa, rOa, my_particles) 05185 ALLOCATE(wcai(SIZE(wc(1)%WannierHamDiag)),stat=ierror) 05186 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05187 ALLOCATE(wcdi(SIZE(wc(1)%WannierHamDiag)),stat=ierror) 05188 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05189 nwca=0 05190 nwcd=0 05191 DO j=1,SIZE(wc(1)%WannierHamDiag) 05192 x=distance(rOd-wc(1)%centres(:,j)) 05193 y=distance(rOa-wc(1)%centres(:,j)) 05194 IF (x<rcut) THEN 05195 nwcd=nwcd+1 05196 wcdi(nwcd)=j 05197 CYCLE 05198 ENDIF 05199 IF (y<rcut) THEN 05200 nwca=nwca+1 05201 wcai(nwca)=j 05202 ENDIF 05203 ENDDO 05204 05205 dmin=distance(rH-wc(1)%centres(:,wcdi(1))) 05206 amin=distance(rH-wc(1)%centres(:,wcai(1))) 05207 idmin=wcdi(1) 05208 iamin=wcai(1) 05209 !dmin constains the smallest numer, amin the next smallest 05210 DO i=2,nwcd 05211 x=distance(rH-wc(1)%centres(:,wcdi(i))) 05212 IF (x<dmin) THEN 05213 dmin=x 05214 idmin=wcdi(i) 05215 ENDIF 05216 ENDDO 05217 DO i=2,nwca 05218 x=distance(rH-wc(1)%centres(:,wcai(i))) 05219 IF (x<amin) THEN 05220 amin=x 05221 iamin=wcai(i) 05222 ENDIF 05223 ENDDO 05224 ! zero=0.0_dp 05225 ! CALL put_derivative(colvar, 1, zero) 05226 ! CALL put_derivative(colvar, 2,zero) 05227 ! CALL put_derivative(colvar, 3, zero) 05228 05229 ! write(6,'(2(i0,1x),4(f16.8,1x))')idmin,iamin,wc(1)%WannierHamDiag(idmin),wc(1)%WannierHamDiag(iamin),dmin,amin 05230 colvar%ss = wc(1)%WannierHamDiag(idmin)-wc(1)%WannierHamDiag(iamin) 05231 DEALLOCATE(wcai,stat=ierror) 05232 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05233 DEALLOCATE(wcdi,stat=ierror) 05234 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05235 05236 CONTAINS 05237 REAL(dp) FUNCTION distance(rij) 05238 REAL(dp), INTENT(in) :: rij(3) 05239 05240 s=MATMUL(cell%h_inv,rij) 05241 s=s-NINT(s) 05242 xv=MATMUL(cell%hmat,s) 05243 distance=SQRT(DOT_PRODUCT(xv,xv)) 05244 END FUNCTION distance 05245 05246 END SUBROUTINE Wc_colvar 05247 05248 ! ***************************************************************************** 05252 SUBROUTINE HBP_colvar(colvar,cell,subsys,particles,qs_env,error) 05253 TYPE(colvar_type), POINTER :: colvar 05254 TYPE(cell_type), POINTER :: cell 05255 TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys 05256 TYPE(particle_type), DIMENSION(:), 05257 OPTIONAL, POINTER :: particles 05258 TYPE(qs_environment_type), POINTER, OPTIONAL :: qs_env ! optional just because I am lazy... but I should get rid of it... 05259 TYPE(cp_error_type), INTENT(inout) :: error 05260 05261 CHARACTER(len=*), PARAMETER :: routineN = 'HBP_colvar', 05262 routineP = moduleN//':'//routineN 05263 05264 INTEGER :: Od, H, Oa 05265 LOGICAL :: failure 05266 REAL(dp) :: rOd(3), rOa(3), rH(3), 05267 x,y,s(3),xv(3),dmin,amin,aux 05268 INTEGER :: idmin, iamin,i,j,ai,il,output_unit 05269 TYPE(particle_list_type), POINTER :: particles_i 05270 TYPE(particle_type), DIMENSION(:), 05271 POINTER :: my_particles 05272 TYPE(wannier_centres_type), 05273 DIMENSION(:),POINTER :: wc 05274 INTEGER, ALLOCATABLE :: wcai(:),wcdi(:) ! contains the indeces of the wannier centres closed to the donor and acceptor 05275 INTEGER :: nwca,nwcd,ierror 05276 REAL(dp) :: rcut 05277 ! real(dp) :: zero(3) 05278 TYPE(cp_logger_type), POINTER :: logger 05279 05280 05281 05282 failure=.FALSE. 05283 NULLIFY(particles_i,wc,logger) 05284 logger => cp_error_get_logger(error) 05285 output_unit= cp_logger_get_default_io_unit(logger) 05286 05287 CPPrecondition(colvar%type_id==HBP_colvar_id,cp_failure_level,routineP,error,failure) 05288 IF (PRESENT(particles)) THEN 05289 my_particles => particles 05290 ELSE 05291 CPPrecondition(PRESENT(subsys),cp_failure_level,routineP,error,failure) 05292 CALL cp_subsys_get(subsys,particles=particles_i,error=error) 05293 my_particles => particles_i%els 05294 END IF 05295 CALL get_qs_env(qs_env,WannierCentres=wc,error=error) 05296 rcut = colvar%HBP%rcut ! distances are in bohr as far as I remember 05297 ALLOCATE(wcai(SIZE(wc(1)%WannierHamDiag)),stat=ierror) 05298 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05299 ALLOCATE(wcdi(SIZE(wc(1)%WannierHamDiag)),stat=ierror) 05300 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05301 colvar%ss=0.0_dp 05302 DO il=1,colvar%HBP%nPoints 05303 Od=colvar%HBP%ids(il,1) 05304 H=colvar%HBP%ids(il,2) 05305 Oa=colvar%HBP%ids(il,3) 05306 CALL get_coordinates(colvar, Od, rOd, my_particles) 05307 CALL get_coordinates(colvar, H, rH, my_particles) 05308 CALL get_coordinates(colvar, Oa, rOa, my_particles) 05309 nwca=0 05310 nwcd=0 05311 DO j=1,SIZE(wc(1)%WannierHamDiag) 05312 x=distance(rOd-wc(1)%centres(:,j)) 05313 y=distance(rOa-wc(1)%centres(:,j)) 05314 IF (x<rcut) THEN 05315 nwcd=nwcd+1 05316 wcdi(nwcd)=j 05317 CYCLE 05318 ENDIF 05319 IF (y<rcut) THEN 05320 nwca=nwca+1 05321 wcai(nwca)=j 05322 ENDIF 05323 ENDDO 05324 05325 dmin=distance(rH-wc(1)%centres(:,wcdi(1))) 05326 amin=distance(rH-wc(1)%centres(:,wcai(1))) 05327 idmin=wcdi(1) 05328 iamin=wcai(1) 05329 !dmin constains the smallest numer, amin the next smallest 05330 DO i=2,nwcd 05331 x=distance(rH-wc(1)%centres(:,wcdi(i))) 05332 IF (x<dmin) THEN 05333 dmin=x 05334 idmin=wcdi(i) 05335 ENDIF 05336 ENDDO 05337 DO i=2,nwca 05338 x=distance(rH-wc(1)%centres(:,wcai(i))) 05339 IF (x<amin) THEN 05340 amin=x 05341 iamin=wcai(i) 05342 ENDIF 05343 ENDDO 05344 ! if (output_unit>0) then 05345 ! write(output_unit,'(a,2(i0,1x),5(f16.8,1x))')"HBP|",idmin,iamin,wc(1)%WannierHamDiag(idmin),wc(1)%WannierHamDiag(iamin),dmin,amin, wc(1)%WannierHamDiag(idmin)-wc(1)%WannierHamDiag(iamin) 05346 ! endif 05347 colvar%HBP%ewc(il) = colvar%HBP%shift + wc(1)%WannierHamDiag(idmin)-wc(1)%WannierHamDiag(iamin) 05348 colvar%ss = colvar%ss + colvar%HBP%shift + wc(1)%WannierHamDiag(idmin)-wc(1)%WannierHamDiag(iamin) 05349 ENDDO 05350 IF (output_unit>0) THEN 05351 DO il=1,colvar%HBP%nPoints 05352 WRITE(output_unit,'(a,1(f16.8,1x))')"HBP| = ",colvar%HBP%ewc(il) 05353 ENDDO 05354 WRITE(output_unit,'(a,1(f16.8,1x))')"HBP|\theta(x) = ",colvar%ss 05355 ENDIF 05356 DEALLOCATE(wcai,stat=ierror) 05357 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05358 DEALLOCATE(wcdi,stat=ierror) 05359 CPPrecondition(ierror==0,cp_failure_level,routineP,error,failure) 05360 05361 CONTAINS 05362 REAL(dp) FUNCTION distance(rij) 05363 REAL(dp), INTENT(in) :: rij(3) 05364 05365 s=MATMUL(cell%h_inv,rij) 05366 s=s-NINT(s) 05367 xv=MATMUL(cell%hmat,s) 05368 distance=SQRT(DOT_PRODUCT(xv,xv)) 05369 END FUNCTION distance 05370 05371 END SUBROUTINE HBP_colvar 05372 05373 05374 END MODULE colvar_methods
1.7.3