|
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 molecule_kind_types 00015 USE atomic_kind_types, ONLY: atomic_kind_type,& 00016 get_atomic_kind 00017 USE basis_set_types, ONLY: get_gto_basis_set,& 00018 gto_basis_set_type 00019 USE colvar_types, ONLY: & 00020 angle_colvar_id, colvar_counters, combine_colvar_id, coord_colvar_id, & 00021 dfunct_colvar_id, dist_colvar_id, gyration_colvar_id, & 00022 hydronium_colvar_id, plane_distance_colvar_id, & 00023 plane_plane_angle_colvar_id, population_colvar_id, qparm_colvar_id, & 00024 reaction_path_colvar_id, rotation_colvar_id, torsion_colvar_id, & 00025 xyz_diag_colvar_id, xyz_outerdiag_colvar_id 00026 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00027 cp_print_key_unit_nr 00028 USE external_potential_types, ONLY: all_potential_type,& 00029 get_potential,& 00030 gth_potential_type 00031 USE f77_blas 00032 USE force_field_types, ONLY: & 00033 bend_kind_type, bond_kind_type, impr_kind_dealloc_ref, impr_kind_type, & 00034 opbend_kind_type, torsion_kind_dealloc_ref, torsion_kind_type, & 00035 ub_kind_dealloc_ref, ub_kind_type 00036 USE input_constants, ONLY: use_perd_x,& 00037 use_perd_xy,& 00038 use_perd_xyz,& 00039 use_perd_xz,& 00040 use_perd_y,& 00041 use_perd_yz,& 00042 use_perd_z 00043 USE input_section_types, ONLY: section_vals_type 00044 USE kinds, ONLY: default_string_length,& 00045 dp 00046 USE shell_potential_types, ONLY: shell_kind_type 00047 USE termination, ONLY: stop_program 00048 USE timings, ONLY: timeset,& 00049 timestop 00050 #include "cp_common_uses.h" 00051 00052 IMPLICIT NONE 00053 00054 PRIVATE 00055 00056 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_kind_types' 00057 00058 ! *** Define the derived structure types *** 00059 00060 ! ***************************************************************************** 00061 TYPE atom_type 00062 TYPE(atomic_kind_type), POINTER :: atomic_kind 00063 INTEGER :: id_name 00064 END TYPE atom_type 00065 00066 ! ***************************************************************************** 00067 TYPE shell_type 00068 INTEGER :: a 00069 CHARACTER(LEN=default_string_length) :: name 00070 TYPE(shell_kind_type), POINTER :: shell_kind 00071 END TYPE shell_type 00072 00073 ! ***************************************************************************** 00074 TYPE bond_type 00075 INTEGER :: a,b 00076 INTEGER :: id_type, itype 00077 TYPE(bond_kind_type), POINTER :: bond_kind 00078 END TYPE bond_type 00079 00080 ! ***************************************************************************** 00081 TYPE bend_type 00082 INTEGER :: a,b,c 00083 INTEGER :: id_type, itype 00084 TYPE(bend_kind_type), POINTER :: bend_kind 00085 END TYPE bend_type 00086 00087 ! ***************************************************************************** 00088 TYPE ub_type 00089 INTEGER :: a,b,c 00090 INTEGER :: id_type, itype 00091 TYPE(ub_kind_type), POINTER :: ub_kind 00092 END TYPE ub_type 00093 00094 ! ***************************************************************************** 00095 TYPE torsion_type 00096 INTEGER :: a,b,c,d 00097 INTEGER :: id_type, itype 00098 TYPE(torsion_kind_type), POINTER :: torsion_kind 00099 END TYPE torsion_type 00100 00101 ! ***************************************************************************** 00102 TYPE impr_type 00103 INTEGER :: a,b,c,d 00104 INTEGER :: id_type, itype 00105 TYPE(impr_kind_type), POINTER :: impr_kind 00106 END TYPE impr_type 00107 00108 ! ***************************************************************************** 00109 TYPE opbend_type 00110 INTEGER :: a,b,c,d 00111 INTEGER :: id_type, itype 00112 TYPE(opbend_kind_type), POINTER :: opbend_kind 00113 END TYPE opbend_type 00114 00115 ! ***************************************************************************** 00116 TYPE restraint_type 00117 LOGICAL :: active 00118 REAL(KIND=dp) :: k0 00119 END TYPE restraint_type 00120 00121 ! ***************************************************************************** 00122 TYPE colvar_constraint_type 00123 INTEGER :: type_id 00124 INTEGER :: inp_seq_num 00125 LOGICAL :: use_points 00126 REAL(KIND = dp) :: expected_value 00127 REAL(KIND = dp) :: expected_value_growth_speed 00128 INTEGER, POINTER, DIMENSION(:) :: i_atoms 00129 TYPE(restraint_type) :: restraint 00130 END TYPE colvar_constraint_type 00131 00132 ! ***************************************************************************** 00133 TYPE g3x3_constraint_type 00134 INTEGER :: a,b,c 00135 REAL(KIND = dp) :: dab,dac,dbc 00136 TYPE(restraint_type) :: restraint 00137 END TYPE g3x3_constraint_type 00138 00139 ! ***************************************************************************** 00140 TYPE g4x6_constraint_type 00141 INTEGER :: a,b,c,d 00142 REAL(KIND = dp) :: dab,dac,dbc,dad,dbd,dcd 00143 TYPE(restraint_type) :: restraint 00144 END TYPE g4x6_constraint_type 00145 00146 ! ***************************************************************************** 00147 TYPE vsite_constraint_type 00148 INTEGER :: a,b,c,d 00149 REAL(KIND = dp) :: wbc, wdc 00150 TYPE(restraint_type) :: restraint 00151 END TYPE vsite_constraint_type 00152 00153 ! ***************************************************************************** 00154 TYPE fixd_constraint_type 00155 TYPE(restraint_type) :: restraint 00156 INTEGER :: fixd, itype 00157 REAL(KIND=dp), DIMENSION(3) :: coord 00158 END TYPE fixd_constraint_type 00159 00160 ! ***************************************************************************** 00161 TYPE local_fixd_constraint_type 00162 INTEGER :: ifixd_index, ikind 00163 END TYPE local_fixd_constraint_type 00164 00165 ! ***************************************************************************** 00166 TYPE molecule_kind_type 00167 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 00168 TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set 00169 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 00170 TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set 00171 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 00172 TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set 00173 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 00174 TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set 00175 TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list 00176 TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set 00177 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 00178 TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set 00179 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 00180 TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list 00181 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list 00182 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list 00183 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list 00184 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 00185 TYPE(shell_type), DIMENSION(:), POINTER :: shell_list 00186 CHARACTER(LEN=default_string_length) :: name 00187 REAL(KIND = dp) :: charge, 00188 mass 00189 INTEGER :: kind_number, 00190 natom, 00191 nbond, 00192 nbend, 00193 nimpr, 00194 nopbend, 00195 ntorsion, 00196 nub, 00197 ng3x3, ng3x3_restraint, 00198 ng4x6, ng4x6_restraint, 00199 nvsite, nvsite_restraint, 00200 nfixd, nfixd_restraint, 00201 nmolecule, nshell 00202 TYPE(colvar_counters) :: ncolv 00203 INTEGER :: nsgf, nelectron 00204 INTEGER, DIMENSION(:), POINTER :: molecule_list 00205 LOGICAL :: molname_generated 00206 END TYPE molecule_kind_type 00207 00208 ! *** Public subroutines *** 00209 PUBLIC :: allocate_molecule_kind_set,& 00210 deallocate_molecule_kind_set,& 00211 get_molecule_kind,& 00212 get_molecule_kind_set,& 00213 num_ao_el_per_molecule,& 00214 set_molecule_kind,& 00215 write_molecule_kind,& 00216 write_molecule_kind_set,& 00217 setup_colvar_counters 00218 00219 ! *** Public data types *** 00220 PUBLIC :: atom_type,& 00221 bend_type,& 00222 bond_type,& 00223 ub_type,& 00224 torsion_type,& 00225 impr_type,& 00226 opbend_type,& 00227 restraint_type,& 00228 colvar_constraint_type,& 00229 g3x3_constraint_type,& 00230 g4x6_constraint_type,& 00231 vsite_constraint_type,& 00232 fixd_constraint_type,& 00233 local_fixd_constraint_type,& 00234 molecule_kind_type,& 00235 shell_type 00236 00237 CONTAINS 00238 00239 ! ***************************************************************************** 00240 SUBROUTINE setup_colvar_counters(colv_list, ncolv,error) 00241 TYPE(colvar_constraint_type), 00242 DIMENSION(:), POINTER :: colv_list 00243 TYPE(colvar_counters) :: ncolv 00244 TYPE(cp_error_type), INTENT(inout) :: error 00245 00246 CHARACTER(len=*), PARAMETER :: routineN = 'setup_colvar_counters', 00247 routineP = moduleN//':'//routineN 00248 00249 INTEGER :: k 00250 LOGICAL :: failure 00251 00252 failure=.FALSE. 00253 ncolv%ndist = 0 00254 ncolv%nangle = 0 00255 ncolv%ndfunct = 0 00256 ncolv%ntorsion = 0 00257 ncolv%ncoord = 0 00258 ncolv%nplane_dist = 0 00259 ncolv%nplane_angle = 0 00260 ncolv%nrot = 0 00261 ncolv%nqparm = 0 00262 ncolv%nxyz_diag = 0 00263 ncolv%nxyz_outerdiag = 0 00264 ncolv%nhydronium = 0 00265 ncolv%nreactionpath = 0 00266 ncolv%ncombinecvs = 0 00267 ncolv%nrestraint = 0 00268 ncolv%npopulation = 0 00269 ncolv%ngyration = 0 00270 00271 IF (ASSOCIATED(colv_list)) THEN 00272 DO k = 1, SIZE(colv_list) 00273 IF (colv_list(k)%restraint%active) ncolv%nrestraint = ncolv%nrestraint + 1 00274 SELECT CASE(colv_list(k)%type_id) 00275 CASE(angle_colvar_id) 00276 ncolv%nangle = ncolv%nangle + 1 00277 CASE(coord_colvar_id) 00278 ncolv%ncoord = ncolv%ncoord + 1 00279 CASE(population_colvar_id) 00280 ncolv%npopulation = ncolv%npopulation + 1 00281 CASE(gyration_colvar_id) 00282 ncolv%ngyration = ncolv%ngyration + 1 00283 CASE(rotation_colvar_id) 00284 ncolv%nrot = ncolv%nrot + 1 00285 CASE(dist_colvar_id) 00286 ncolv%ndist = ncolv%ndist + 1 00287 CASE(dfunct_colvar_id) 00288 ncolv%ndfunct = ncolv%ndfunct + 1 00289 CASE(plane_distance_colvar_id) 00290 ncolv%nplane_dist = ncolv%nplane_dist + 1 00291 CASE(plane_plane_angle_colvar_id) 00292 ncolv%nplane_angle = ncolv%nplane_angle + 1 00293 CASE(torsion_colvar_id) 00294 ncolv%ntorsion = ncolv%ntorsion + 1 00295 CASE(qparm_colvar_id) 00296 ncolv%nqparm = ncolv%nqparm + 1 00297 CASE(xyz_diag_colvar_id) 00298 ncolv%nxyz_diag = ncolv%nxyz_diag + 1 00299 CASE(xyz_outerdiag_colvar_id) 00300 ncolv%nxyz_outerdiag = ncolv%nxyz_outerdiag + 1 00301 CASE(hydronium_colvar_id) 00302 ncolv%nhydronium = ncolv%nhydronium + 1 00303 CASE(reaction_path_colvar_id) 00304 ncolv%nreactionpath = ncolv%nreactionpath + 1 00305 CASE(combine_colvar_id) 00306 ncolv%ncombinecvs = ncolv%ncombinecvs + 1 00307 CASE DEFAULT 00308 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00309 END SELECT 00310 END DO 00311 END IF 00312 ncolv%ntot = ncolv%ndist +& 00313 ncolv%nangle +& 00314 ncolv%ntorsion +& 00315 ncolv%ncoord +& 00316 ncolv%nplane_dist +& 00317 ncolv%nplane_angle +& 00318 ncolv%ndfunct +& 00319 ncolv%nrot +& 00320 ncolv%nqparm +& 00321 ncolv%nxyz_diag +& 00322 ncolv%nxyz_outerdiag +& 00323 ncolv%nhydronium +& 00324 ncolv%nreactionpath +& 00325 ncolv%ncombinecvs +& 00326 ncolv%npopulation +& 00327 ncolv%ngyration 00328 00329 00330 END SUBROUTINE setup_colvar_counters 00331 00332 ! ***************************************************************************** 00338 SUBROUTINE allocate_molecule_kind_set(molecule_kind_set,nmolecule_kind,error) 00339 TYPE(molecule_kind_type), DIMENSION(:), 00340 POINTER :: molecule_kind_set 00341 INTEGER, INTENT(IN) :: nmolecule_kind 00342 TYPE(cp_error_type), INTENT(inout) :: error 00343 00344 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_molecule_kind_set', 00345 routineP = moduleN//':'//routineN 00346 00347 INTEGER :: imolecule_kind, stat 00348 LOGICAL :: failure 00349 00350 failure = .FALSE. 00351 IF (ASSOCIATED(molecule_kind_set)) THEN 00352 CALL deallocate_molecule_kind_set(molecule_kind_set,error=error) 00353 END IF 00354 00355 ALLOCATE (molecule_kind_set(nmolecule_kind),STAT=stat) 00356 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00357 00358 DO imolecule_kind=1,nmolecule_kind 00359 NULLIFY (molecule_kind_set(imolecule_kind)%atom_list) 00360 NULLIFY (molecule_kind_set(imolecule_kind)%bond_list) 00361 NULLIFY (molecule_kind_set(imolecule_kind)%bend_list) 00362 NULLIFY (molecule_kind_set(imolecule_kind)%colv_list) 00363 NULLIFY (molecule_kind_set(imolecule_kind)%ub_list) 00364 NULLIFY (molecule_kind_set(imolecule_kind)%ub_kind_set) 00365 NULLIFY (molecule_kind_set(imolecule_kind)%impr_kind_set) 00366 NULLIFY (molecule_kind_set(imolecule_kind)%impr_list) 00367 NULLIFY (molecule_kind_set(imolecule_kind)%opbend_kind_set) 00368 NULLIFY (molecule_kind_set(imolecule_kind)%opbend_list) 00369 NULLIFY (molecule_kind_set(imolecule_kind)%g3x3_list) 00370 NULLIFY (molecule_kind_set(imolecule_kind)%g4x6_list) 00371 NULLIFY (molecule_kind_set(imolecule_kind)%vsite_list) 00372 NULLIFY (molecule_kind_set(imolecule_kind)%fixd_list) 00373 NULLIFY (molecule_kind_set(imolecule_kind)%shell_list) 00374 NULLIFY (molecule_kind_set(imolecule_kind)%torsion_list) 00375 NULLIFY (molecule_kind_set(imolecule_kind)%bond_kind_set) 00376 NULLIFY (molecule_kind_set(imolecule_kind)%bend_kind_set) 00377 NULLIFY (molecule_kind_set(imolecule_kind)%torsion_kind_set) 00378 molecule_kind_set(imolecule_kind)%charge = 0.0_dp 00379 molecule_kind_set(imolecule_kind)%mass = 0.0_dp 00380 molecule_kind_set(imolecule_kind)%name = "" 00381 molecule_kind_set(imolecule_kind)%molname_generated = .FALSE. 00382 molecule_kind_set(imolecule_kind)%kind_number = imolecule_kind 00383 molecule_kind_set(imolecule_kind)%natom = 0 00384 molecule_kind_set(imolecule_kind)%nbend = 0 00385 molecule_kind_set(imolecule_kind)%nbond = 0 00386 molecule_kind_set(imolecule_kind)%nimpr = 0 00387 molecule_kind_set(imolecule_kind)%nopbend = 0 00388 molecule_kind_set(imolecule_kind)%nub = 0 00389 CALL setup_colvar_counters(molecule_kind_set(imolecule_kind)%colv_list,& 00390 molecule_kind_set(imolecule_kind)%ncolv,error) 00391 molecule_kind_set(imolecule_kind)%ng3x3 = 0 00392 molecule_kind_set(imolecule_kind)%ng4x6 = 0 00393 molecule_kind_set(imolecule_kind)%nvsite = 0 00394 molecule_kind_set(imolecule_kind)%nfixd = 0 00395 molecule_kind_set(imolecule_kind)%ng3x3_restraint = 0 00396 molecule_kind_set(imolecule_kind)%ng4x6_restraint = 0 00397 molecule_kind_set(imolecule_kind)%nvsite_restraint = 0 00398 molecule_kind_set(imolecule_kind)%nfixd_restraint = 0 00399 molecule_kind_set(imolecule_kind)%nmolecule = 0 00400 molecule_kind_set(imolecule_kind)%ntorsion = 0 00401 molecule_kind_set(imolecule_kind)%nshell = 0 00402 NULLIFY (molecule_kind_set(imolecule_kind)%molecule_list) 00403 END DO 00404 00405 END SUBROUTINE allocate_molecule_kind_set 00406 00407 ! ***************************************************************************** 00413 SUBROUTINE deallocate_molecule_kind_set(molecule_kind_set,error) 00414 00415 TYPE(molecule_kind_type), DIMENSION(:), 00416 POINTER :: molecule_kind_set 00417 TYPE(cp_error_type), INTENT(inout) :: error 00418 00419 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_molecule_kind_set', 00420 routineP = moduleN//':'//routineN 00421 00422 INTEGER :: i, imolecule_kind, j, 00423 nmolecule_kind, stat 00424 LOGICAL :: failure 00425 00426 failure = .FALSE. 00427 IF (ASSOCIATED(molecule_kind_set)) THEN 00428 00429 nmolecule_kind = SIZE(molecule_kind_set) 00430 00431 DO imolecule_kind=1,nmolecule_kind 00432 00433 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%atom_list)) THEN 00434 DEALLOCATE (molecule_kind_set(imolecule_kind)%atom_list,STAT=stat) 00435 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00436 END IF 00437 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bend_kind_set)) THEN 00438 DEALLOCATE (molecule_kind_set(imolecule_kind)%bend_kind_set,STAT=stat) 00439 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00440 END IF 00441 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bend_list)) THEN 00442 DEALLOCATE (molecule_kind_set(imolecule_kind)%bend_list,STAT=stat) 00443 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00444 END IF 00445 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%ub_list)) THEN 00446 DEALLOCATE (molecule_kind_set(imolecule_kind)%ub_list,STAT=stat) 00447 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00448 END IF 00449 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%ub_kind_set)) THEN 00450 CALL ub_kind_dealloc_ref(molecule_kind_set(imolecule_kind)%ub_kind_set,error=error) 00451 END IF 00452 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%impr_list)) THEN 00453 DEALLOCATE (molecule_kind_set(imolecule_kind)%impr_list,STAT=stat) 00454 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00455 END IF 00456 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%impr_kind_set)) THEN 00457 DO i = 1, SIZE(molecule_kind_set(imolecule_kind)%impr_kind_set) 00458 CALL impr_kind_dealloc_ref(molecule_kind_set(imolecule_kind)%impr_kind_set(i),error=error) 00459 END DO 00460 DEALLOCATE (molecule_kind_set(imolecule_kind)%impr_kind_set,STAT=stat) 00461 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00462 END IF 00463 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_list)) THEN 00464 DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_list,STAT=stat) 00465 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00466 END IF 00467 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_kind_set)) THEN 00468 DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_kind_set,STAT=stat) 00469 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00470 END IF 00471 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bond_kind_set)) THEN 00472 DEALLOCATE (molecule_kind_set(imolecule_kind)%bond_kind_set,STAT=stat) 00473 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00474 END IF 00475 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bond_list)) THEN 00476 DEALLOCATE (molecule_kind_set(imolecule_kind)%bond_list,STAT=stat) 00477 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00478 END IF 00479 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%colv_list)) THEN 00480 DO j = 1, SIZE(molecule_kind_set(imolecule_kind)%colv_list) 00481 DEALLOCATE (molecule_kind_set(imolecule_kind)%colv_list(j)%i_atoms,STAT=stat) 00482 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00483 END DO 00484 DEALLOCATE (molecule_kind_set(imolecule_kind)%colv_list,STAT=stat) 00485 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00486 END IF 00487 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%g3x3_list)) THEN 00488 DEALLOCATE (molecule_kind_set(imolecule_kind)%g3x3_list,STAT=stat) 00489 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00490 END IF 00491 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%g4x6_list)) THEN 00492 DEALLOCATE (molecule_kind_set(imolecule_kind)%g4x6_list,STAT=stat) 00493 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00494 END IF 00495 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%vsite_list)) THEN 00496 DEALLOCATE (molecule_kind_set(imolecule_kind)%vsite_list,STAT=stat) 00497 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00498 END IF 00499 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%fixd_list)) THEN 00500 DEALLOCATE (molecule_kind_set(imolecule_kind)%fixd_list,STAT=stat) 00501 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00502 END IF 00503 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%torsion_kind_set)) THEN 00504 DO i = 1, SIZE(molecule_kind_set(imolecule_kind)%torsion_kind_set) 00505 CALL torsion_kind_dealloc_ref(molecule_kind_set(imolecule_kind)%torsion_kind_set(i),error=error) 00506 END DO 00507 DEALLOCATE (molecule_kind_set(imolecule_kind)%torsion_kind_set,STAT=stat) 00508 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00509 END IF 00510 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%shell_list)) THEN 00511 DEALLOCATE (molecule_kind_set(imolecule_kind)%shell_list,STAT=stat) 00512 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00513 END IF 00514 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%torsion_list)) THEN 00515 DEALLOCATE (molecule_kind_set(imolecule_kind)%torsion_list,STAT=stat) 00516 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00517 END IF 00518 IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%molecule_list)) THEN 00519 DEALLOCATE (molecule_kind_set(imolecule_kind)%molecule_list,STAT=stat) 00520 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00521 ENDIF 00522 END DO 00523 00524 DEALLOCATE (molecule_kind_set,STAT=stat) 00525 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00526 ELSE 00527 CALL stop_program(routineN,moduleN,__LINE__,& 00528 "The pointer molecule_kind_set is not associated and cannot be deallocated") 00529 END IF 00530 00531 END SUBROUTINE deallocate_molecule_kind_set 00532 00533 ! ***************************************************************************** 00539 SUBROUTINE get_molecule_kind(molecule_kind,atom_list,bond_list,bend_list,& 00540 ub_list,impr_list,opbend_list,colv_list,fixd_list,& 00541 g3x3_list,g4x6_list,vsite_list,torsion_list,shell_list,& 00542 name,mass,charge,kind_number,natom,nbend,nbond,nub,& 00543 nimpr,nopbend,nconstraint,nconstraint_fixd,nfixd,ncolv,ng3x3,ng4x6,& 00544 nvsite,nfixd_restraint,ng3x3_restraint, ng4x6_restraint,& 00545 nvsite_restraint,nrestraints,nmolecule,nsgf,nshell,ntorsion,& 00546 molecule_list,nelectron, bond_kind_set,bend_kind_set,& 00547 ub_kind_set,impr_kind_set,opbend_kind_set,torsion_kind_set,& 00548 molname_generated) 00549 00550 TYPE(molecule_kind_type), POINTER :: molecule_kind 00551 TYPE(atom_type), DIMENSION(:), 00552 OPTIONAL, POINTER :: atom_list 00553 TYPE(bond_type), DIMENSION(:), 00554 OPTIONAL, POINTER :: bond_list 00555 TYPE(bend_type), DIMENSION(:), 00556 OPTIONAL, POINTER :: bend_list 00557 TYPE(ub_type), DIMENSION(:), OPTIONAL, 00558 POINTER :: ub_list 00559 TYPE(impr_type), DIMENSION(:), 00560 OPTIONAL, POINTER :: impr_list 00561 TYPE(opbend_type), DIMENSION(:), 00562 OPTIONAL, POINTER :: opbend_list 00563 TYPE(colvar_constraint_type), 00564 DIMENSION(:), OPTIONAL, POINTER :: colv_list 00565 TYPE(fixd_constraint_type), 00566 DIMENSION(:), OPTIONAL, POINTER :: fixd_list 00567 TYPE(g3x3_constraint_type), 00568 DIMENSION(:), OPTIONAL, POINTER :: g3x3_list 00569 TYPE(g4x6_constraint_type), 00570 DIMENSION(:), OPTIONAL, POINTER :: g4x6_list 00571 TYPE(vsite_constraint_type), 00572 DIMENSION(:), OPTIONAL, POINTER :: vsite_list 00573 TYPE(torsion_type), DIMENSION(:), 00574 OPTIONAL, POINTER :: torsion_list 00575 TYPE(shell_type), DIMENSION(:), 00576 OPTIONAL, POINTER :: shell_list 00577 CHARACTER(LEN=default_string_length), 00578 INTENT(OUT), OPTIONAL :: name 00579 REAL(KIND=dp), OPTIONAL :: mass, charge 00580 INTEGER, INTENT(OUT), OPTIONAL :: kind_number, natom, nbend, 00581 nbond, nub, nimpr, nopbend, 00582 nconstraint, 00583 nconstraint_fixd, nfixd 00584 TYPE(colvar_counters), INTENT(out), 00585 OPTIONAL :: ncolv 00586 INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nvsite, nfixd_restraint, 00587 ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, 00588 nmolecule, nsgf, nshell, ntorsion 00589 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list 00590 INTEGER, INTENT(OUT), OPTIONAL :: nelectron 00591 TYPE(bond_kind_type), DIMENSION(:), 00592 OPTIONAL, POINTER :: bond_kind_set 00593 TYPE(bend_kind_type), DIMENSION(:), 00594 OPTIONAL, POINTER :: bend_kind_set 00595 TYPE(ub_kind_type), DIMENSION(:), 00596 OPTIONAL, POINTER :: ub_kind_set 00597 TYPE(impr_kind_type), DIMENSION(:), 00598 OPTIONAL, POINTER :: impr_kind_set 00599 TYPE(opbend_kind_type), DIMENSION(:), 00600 OPTIONAL, POINTER :: opbend_kind_set 00601 TYPE(torsion_kind_type), DIMENSION(:), 00602 OPTIONAL, POINTER :: torsion_kind_set 00603 LOGICAL, INTENT(OUT), OPTIONAL :: molname_generated 00604 00605 CHARACTER(len=*), PARAMETER :: routineN = 'get_molecule_kind', 00606 routineP = moduleN//':'//routineN 00607 00608 INTEGER :: i 00609 00610 IF (ASSOCIATED(molecule_kind)) THEN 00611 00612 IF (PRESENT(atom_list)) atom_list => molecule_kind%atom_list 00613 IF (PRESENT(bend_list)) bend_list => molecule_kind%bend_list 00614 IF (PRESENT(bond_list)) bond_list => molecule_kind%bond_list 00615 IF (PRESENT(impr_list)) impr_list => molecule_kind%impr_list 00616 IF (PRESENT(opbend_list)) opbend_list => molecule_kind%opbend_list 00617 IF (PRESENT(ub_list)) ub_list => molecule_kind%ub_list 00618 IF (PRESENT(bond_kind_set)) bond_kind_set => molecule_kind%bond_kind_set 00619 IF (PRESENT(bend_kind_set)) bend_kind_set => molecule_kind%bend_kind_set 00620 IF (PRESENT(ub_kind_set)) ub_kind_set => molecule_kind%ub_kind_set 00621 IF (PRESENT(impr_kind_set)) impr_kind_set => molecule_kind%impr_kind_set 00622 IF (PRESENT(opbend_kind_set)) opbend_kind_set => molecule_kind%opbend_kind_set 00623 IF (PRESENT(torsion_kind_set)) torsion_kind_set => molecule_kind%torsion_kind_set 00624 IF (PRESENT(colv_list)) colv_list => molecule_kind%colv_list 00625 IF (PRESENT(g3x3_list)) g3x3_list => molecule_kind%g3x3_list 00626 IF (PRESENT(g4x6_list)) g4x6_list => molecule_kind%g4x6_list 00627 IF (PRESENT(vsite_list)) vsite_list => molecule_kind%vsite_list 00628 IF (PRESENT(fixd_list)) fixd_list => molecule_kind%fixd_list 00629 IF (PRESENT(torsion_list)) torsion_list => molecule_kind%torsion_list 00630 IF (PRESENT(shell_list)) shell_list => molecule_kind%shell_list 00631 IF (PRESENT(name)) name = molecule_kind%name 00632 IF (PRESENT(molname_generated)) molname_generated = molecule_kind%molname_generated 00633 IF (PRESENT(mass)) mass = molecule_kind%mass 00634 IF (PRESENT(charge)) charge = molecule_kind%charge 00635 IF (PRESENT(kind_number)) kind_number = molecule_kind%kind_number 00636 IF (PRESENT(natom)) natom = molecule_kind%natom 00637 IF (PRESENT(nbend)) nbend = molecule_kind%nbend 00638 IF (PRESENT(nbond)) nbond = molecule_kind%nbond 00639 IF (PRESENT(nub)) nub = molecule_kind%nub 00640 IF (PRESENT(nimpr)) nimpr = molecule_kind%nimpr 00641 IF (PRESENT(nopbend)) nopbend = molecule_kind%nopbend 00642 IF (PRESENT(nconstraint)) nconstraint = (molecule_kind%ncolv%ntot - molecule_kind%ncolv%nrestraint) +& 00643 3*(molecule_kind%ng3x3-molecule_kind%ng3x3_restraint) +& 00644 6*(molecule_kind%ng4x6-molecule_kind%ng4x6_restraint)+& 00645 3*(molecule_kind%nvsite-molecule_kind%nvsite_restraint) 00646 IF (PRESENT(ncolv)) ncolv = molecule_kind%ncolv 00647 IF (PRESENT(ng3x3)) ng3x3= molecule_kind%ng3x3 00648 IF (PRESENT(ng4x6)) ng4x6= molecule_kind%ng4x6 00649 IF (PRESENT(nvsite)) nvsite= molecule_kind%nvsite 00650 ! Number of atoms that have one or more components fixed 00651 IF (PRESENT(nfixd)) nfixd= molecule_kind%nfixd 00652 ! Number of degrees of freedom fixed 00653 IF (PRESENT(nconstraint_fixd)) THEN 00654 nconstraint_fixd = 0 00655 IF (molecule_kind%nfixd/=0) THEN 00656 DO i = 1, SIZE(molecule_kind%fixd_list) 00657 IF (molecule_kind%fixd_list(i)%restraint%active) CYCLE 00658 SELECT CASE(molecule_kind%fixd_list(i)%itype) 00659 CASE(use_perd_x, use_perd_y, use_perd_z) 00660 nconstraint_fixd=nconstraint_fixd+1 00661 CASE(use_perd_xy, use_perd_xz, use_perd_yz) 00662 nconstraint_fixd=nconstraint_fixd+2 00663 CASE(use_perd_xyz) 00664 nconstraint_fixd=nconstraint_fixd+3 00665 END SELECT 00666 END DO 00667 END IF 00668 END IF 00669 IF (PRESENT(ng3x3_restraint)) ng3x3_restraint= molecule_kind%ng3x3_restraint 00670 IF (PRESENT(ng4x6_restraint)) ng4x6_restraint= molecule_kind%ng4x6_restraint 00671 IF (PRESENT(nvsite_restraint)) nvsite_restraint= molecule_kind%nvsite_restraint 00672 IF (PRESENT(nfixd_restraint)) nfixd_restraint= molecule_kind%nfixd_restraint 00673 IF (PRESENT(nrestraints)) nrestraints = molecule_kind%ncolv%nrestraint +& 00674 molecule_kind%ng3x3_restraint +& 00675 molecule_kind%ng4x6_restraint+& 00676 molecule_kind%nvsite_restraint 00677 IF (PRESENT(nmolecule)) nmolecule = molecule_kind%nmolecule 00678 IF (PRESENT(nshell)) nshell= molecule_kind%nshell 00679 IF (PRESENT(ntorsion)) ntorsion= molecule_kind%ntorsion 00680 IF (PRESENT(nsgf)) nsgf = molecule_kind%nsgf 00681 IF (PRESENT(nelectron)) nelectron = molecule_kind%nelectron 00682 IF (PRESENT(molecule_list)) molecule_list => molecule_kind%molecule_list 00683 00684 ELSE 00685 00686 CALL stop_program(routineN,moduleN,__LINE__,& 00687 "The pointer molecule_kind is not associated") 00688 00689 END IF 00690 00691 END SUBROUTINE get_molecule_kind 00692 00693 ! ***************************************************************************** 00699 SUBROUTINE get_molecule_kind_set(molecule_kind_set,maxatom,natom,& 00700 nbond,nbend,nub,ntorsion,nimpr,nopbend,& 00701 nconstraint,nconstraint_fixd,nmolecule,& 00702 nrestraints) 00703 00704 TYPE(molecule_kind_type), DIMENSION(:), 00705 POINTER :: molecule_kind_set 00706 INTEGER, INTENT(OUT), OPTIONAL :: maxatom, natom, nbond, nbend, nub, 00707 ntorsion, nimpr, nopbend, nconstraint, nconstraint_fixd, nmolecule, 00708 nrestraints 00709 00710 CHARACTER(len=*), PARAMETER :: routineN = 'get_molecule_kind_set', 00711 routineP = moduleN//':'//routineN 00712 00713 INTEGER :: ibend, ibond, iimpr, imolecule_kind, iopbend, itorsion, iub, 00714 na, nc, nc_fixd, nfixd_restraint, nm, nmolecule_kind, nrestraints_tot 00715 TYPE(molecule_kind_type), POINTER :: molecule_kind 00716 00717 IF (ASSOCIATED(molecule_kind_set)) THEN 00718 00719 IF (PRESENT(maxatom)) maxatom = 0 00720 IF (PRESENT(natom)) natom = 0 00721 IF (PRESENT(nbond)) nbond = 0 00722 IF (PRESENT(nbend)) nbend = 0 00723 IF (PRESENT(nub)) nub = 0 00724 IF (PRESENT(ntorsion)) ntorsion = 0 00725 IF (PRESENT(nimpr)) nimpr = 0 00726 IF (PRESENT(nopbend)) nopbend = 0 00727 IF (PRESENT(nconstraint)) nconstraint = 0 00728 IF (PRESENT(nconstraint_fixd)) nconstraint_fixd = 0 00729 IF (PRESENT(nrestraints)) nrestraints = 0 00730 IF (PRESENT(nmolecule)) nmolecule = 0 00731 00732 nmolecule_kind = SIZE(molecule_kind_set) 00733 00734 DO imolecule_kind=1,nmolecule_kind 00735 00736 molecule_kind => molecule_kind_set(imolecule_kind) 00737 00738 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00739 natom=na,& 00740 nbond=ibond,& 00741 nbend=ibend,& 00742 nub=iub,& 00743 ntorsion=itorsion,& 00744 nimpr=iimpr,& 00745 nopbend=iopbend,& 00746 nconstraint=nc,& 00747 nconstraint_fixd=nc_fixd,& 00748 nfixd_restraint=nfixd_restraint,& 00749 nrestraints=nrestraints_tot,& 00750 nmolecule=nm) 00751 IF (PRESENT(maxatom)) maxatom = MAX(maxatom,na) 00752 IF (PRESENT(natom)) natom = natom + na*nm 00753 IF (PRESENT(nbond)) nbond = nbond + ibond*nm 00754 IF (PRESENT(nbend)) nbend = nbend + ibend*nm 00755 IF (PRESENT(nub)) nub = nub + iub*nm 00756 IF (PRESENT(ntorsion)) ntorsion = ntorsion + itorsion*nm 00757 IF (PRESENT(nimpr)) nimpr = nimpr + iimpr*nm 00758 IF (PRESENT(nopbend)) nopbend = nopbend + iopbend*nm 00759 IF (PRESENT(nconstraint)) nconstraint = nconstraint + nc*nm + nc_fixd 00760 IF (PRESENT(nconstraint_fixd)) nconstraint_fixd = nconstraint_fixd + nc_fixd 00761 IF (PRESENT(nmolecule)) nmolecule = nmolecule + nm 00762 IF (PRESENT(nrestraints)) nrestraints = nrestraints + nm*nrestraints_tot + nfixd_restraint 00763 00764 END DO 00765 00766 ELSE 00767 00768 CALL stop_program(routineN,moduleN,__LINE__,& 00769 "The pointer molecule_kind_set is not associated") 00770 00771 END IF 00772 00773 END SUBROUTINE get_molecule_kind_set 00774 00775 ! ***************************************************************************** 00781 SUBROUTINE num_ao_el_per_molecule(molecule_kind_set,charge_x_mol) 00782 00783 TYPE(molecule_kind_type), DIMENSION(:), 00784 POINTER :: molecule_kind_set 00785 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: charge_x_mol 00786 00787 CHARACTER(len=*), PARAMETER :: routineN = 'num_ao_el_per_molecule', 00788 routineP = moduleN//':'//routineN 00789 00790 INTEGER :: iatom, imol, my_charge, n_ao, 00791 natom, nelectron, nmol_kind, 00792 nsgf 00793 REAL(KIND=dp) :: zeff, zeff_correction 00794 TYPE(all_potential_type), POINTER :: all_potential 00795 TYPE(atomic_kind_type), POINTER :: atomic_kind 00796 TYPE(gth_potential_type), POINTER :: gth_potential 00797 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00798 TYPE(molecule_kind_type), POINTER :: molecule_kind 00799 00800 IF (ASSOCIATED(molecule_kind_set)) THEN 00801 00802 nmol_kind = SIZE(molecule_kind_set,1) 00803 natom = 0 00804 my_charge = 0 00805 ! *** Initialize the molecule kind data structure *** 00806 00807 DO imol = 1,nmol_kind 00808 molecule_kind => molecule_kind_set(imol) 00809 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00810 natom=natom) 00811 nelectron = 0 00812 n_ao = 0 00813 IF(PRESENT(charge_x_mol)) THEN 00814 my_charge = charge_x_mol(imol) 00815 END IF 00816 DO iatom=1,natom 00817 atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind 00818 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00819 all_potential=all_potential,& 00820 gth_potential=gth_potential,& 00821 orb_basis_set=orb_basis_set) 00822 IF (ASSOCIATED(all_potential)) THEN 00823 CALL get_potential(potential=all_potential,zeff=zeff,& 00824 zeff_correction=zeff_correction) 00825 ELSE IF (ASSOCIATED(gth_potential)) THEN 00826 CALL get_potential(potential=gth_potential,zeff=zeff,& 00827 zeff_correction=zeff_correction) 00828 ELSE 00829 zeff = 0.0_dp 00830 zeff_correction = 0.0_dp 00831 END IF 00832 nelectron = nelectron + NINT(zeff-zeff_correction) 00833 00834 IF (ASSOCIATED(orb_basis_set)) THEN 00835 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,nsgf=nsgf) 00836 ELSE 00837 nsgf = 0 00838 ENDIF 00839 n_ao = n_ao + nsgf 00840 END DO ! iatom 00841 nelectron = nelectron + my_charge 00842 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00843 nelectron=nelectron, nsgf=n_ao) 00844 END DO ! imol 00845 END IF 00846 00847 END SUBROUTINE num_ao_el_per_molecule 00848 00849 ! ***************************************************************************** 00855 SUBROUTINE set_molecule_kind(molecule_kind,name,mass,charge,kind_number,& 00856 molecule_list,atom_list,nbond,bond_list,& 00857 nbend,bend_list,nub,ub_list,nimpr,impr_list,& 00858 nopbend, opbend_list, ntorsion,& 00859 torsion_list,fixd_list,ncolv,colv_list,ng3x3,& 00860 g3x3_list,ng4x6,nfixd,g4x6_list,nvsite,& 00861 vsite_list,ng3x3_restraint,ng4x6_restraint,& 00862 nfixd_restraint,nshell, shell_list,& 00863 nvsite_restraint,bond_kind_set,bend_kind_set,& 00864 ub_kind_set,torsion_kind_set,impr_kind_set,& 00865 opbend_kind_set, nelectron,nsgf,& 00866 molname_generated) 00867 00868 TYPE(molecule_kind_type), POINTER :: molecule_kind 00869 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name 00870 REAL(KIND=dp), OPTIONAL :: mass, charge 00871 INTEGER, INTENT(IN), OPTIONAL :: kind_number 00872 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list 00873 TYPE(atom_type), DIMENSION(:), 00874 OPTIONAL, POINTER :: atom_list 00875 INTEGER, INTENT(IN), OPTIONAL :: nbond 00876 TYPE(bond_type), DIMENSION(:), 00877 OPTIONAL, POINTER :: bond_list 00878 INTEGER, INTENT(IN), OPTIONAL :: nbend 00879 TYPE(bend_type), DIMENSION(:), 00880 OPTIONAL, POINTER :: bend_list 00881 INTEGER, INTENT(IN), OPTIONAL :: nub 00882 TYPE(ub_type), DIMENSION(:), OPTIONAL, 00883 POINTER :: ub_list 00884 INTEGER, INTENT(IN), OPTIONAL :: nimpr 00885 TYPE(impr_type), DIMENSION(:), 00886 OPTIONAL, POINTER :: impr_list 00887 INTEGER, INTENT(IN), OPTIONAL :: nopbend 00888 TYPE(opbend_type), DIMENSION(:), 00889 OPTIONAL, POINTER :: opbend_list 00890 INTEGER, INTENT(IN), OPTIONAL :: ntorsion 00891 TYPE(torsion_type), DIMENSION(:), 00892 OPTIONAL, POINTER :: torsion_list 00893 TYPE(fixd_constraint_type), 00894 DIMENSION(:), OPTIONAL, POINTER :: fixd_list 00895 TYPE(colvar_counters), INTENT(IN), 00896 OPTIONAL :: ncolv 00897 TYPE(colvar_constraint_type), 00898 DIMENSION(:), OPTIONAL, POINTER :: colv_list 00899 INTEGER, INTENT(IN), OPTIONAL :: ng3x3 00900 TYPE(g3x3_constraint_type), 00901 DIMENSION(:), OPTIONAL, POINTER :: g3x3_list 00902 INTEGER, INTENT(IN), OPTIONAL :: ng4x6, nfixd 00903 TYPE(g4x6_constraint_type), 00904 DIMENSION(:), OPTIONAL, POINTER :: g4x6_list 00905 INTEGER, INTENT(IN), OPTIONAL :: nvsite 00906 TYPE(vsite_constraint_type), 00907 DIMENSION(:), OPTIONAL, POINTER :: vsite_list 00908 INTEGER, INTENT(IN), OPTIONAL :: ng3x3_restraint, 00909 ng4x6_restraint, 00910 nfixd_restraint, nshell 00911 TYPE(shell_type), DIMENSION(:), 00912 OPTIONAL, POINTER :: shell_list 00913 INTEGER, INTENT(IN), OPTIONAL :: nvsite_restraint 00914 TYPE(bond_kind_type), DIMENSION(:), 00915 OPTIONAL, POINTER :: bond_kind_set 00916 TYPE(bend_kind_type), DIMENSION(:), 00917 OPTIONAL, POINTER :: bend_kind_set 00918 TYPE(ub_kind_type), DIMENSION(:), 00919 OPTIONAL, POINTER :: ub_kind_set 00920 TYPE(torsion_kind_type), DIMENSION(:), 00921 OPTIONAL, POINTER :: torsion_kind_set 00922 TYPE(impr_kind_type), DIMENSION(:), 00923 OPTIONAL, POINTER :: impr_kind_set 00924 TYPE(opbend_kind_type), DIMENSION(:), 00925 OPTIONAL, POINTER :: opbend_kind_set 00926 INTEGER, INTENT(IN), OPTIONAL :: nelectron, nsgf 00927 LOGICAL, INTENT(IN), OPTIONAL :: molname_generated 00928 00929 CHARACTER(len=*), PARAMETER :: routineN = 'set_molecule_kind', 00930 routineP = moduleN//':'//routineN 00931 00932 INTEGER :: n 00933 00934 IF (ASSOCIATED(molecule_kind)) THEN 00935 00936 IF (PRESENT(atom_list)) THEN 00937 n = SIZE(atom_list) 00938 molecule_kind%natom = n 00939 molecule_kind%atom_list => atom_list 00940 END IF 00941 IF (PRESENT(molname_generated)) molecule_kind%molname_generated = molname_generated 00942 IF (PRESENT(name)) molecule_kind%name = name 00943 IF (PRESENT(mass)) molecule_kind%mass = mass 00944 IF (PRESENT(charge)) molecule_kind%charge = charge 00945 IF (PRESENT(kind_number)) molecule_kind%kind_number = kind_number 00946 IF (PRESENT(nbond)) molecule_kind%nbond = nbond 00947 IF (PRESENT(bond_list)) molecule_kind%bond_list => bond_list 00948 IF (PRESENT(nbend)) molecule_kind%nbend = nbend 00949 IF (PRESENT(nelectron)) molecule_kind%nelectron = nelectron 00950 IF (PRESENT(nsgf)) molecule_kind%nsgf = nsgf 00951 IF (PRESENT(bend_list)) molecule_kind%bend_list => bend_list 00952 IF (PRESENT(nub)) molecule_kind%nub = nub 00953 IF (PRESENT(ub_list)) molecule_kind%ub_list => ub_list 00954 IF (PRESENT(ntorsion)) molecule_kind%ntorsion = ntorsion 00955 IF (PRESENT(torsion_list)) molecule_kind%torsion_list => torsion_list 00956 IF (PRESENT(nimpr)) molecule_kind%nimpr = nimpr 00957 IF (PRESENT(impr_list)) molecule_kind%impr_list => impr_list 00958 IF (PRESENT(nopbend)) molecule_kind%nopbend = nopbend 00959 IF (PRESENT(opbend_list)) molecule_kind%opbend_list => opbend_list 00960 IF (PRESENT(ncolv)) molecule_kind%ncolv = ncolv 00961 IF (PRESENT(colv_list)) molecule_kind%colv_list => colv_list 00962 IF (PRESENT(ng3x3)) molecule_kind%ng3x3 = ng3x3 00963 IF (PRESENT(g3x3_list)) molecule_kind%g3x3_list => g3x3_list 00964 IF (PRESENT(ng4x6)) molecule_kind%ng4x6 = ng4x6 00965 IF (PRESENT(nvsite)) molecule_kind%nvsite = nvsite 00966 IF (PRESENT(nfixd)) molecule_kind%nfixd = nfixd 00967 IF (PRESENT(nfixd_restraint)) molecule_kind%nfixd_restraint = nfixd_restraint 00968 IF (PRESENT(ng3x3_restraint)) molecule_kind%ng3x3_restraint = ng3x3_restraint 00969 IF (PRESENT(ng4x6_restraint)) molecule_kind%ng4x6_restraint = ng4x6_restraint 00970 IF (PRESENT(nvsite_restraint)) molecule_kind%nvsite_restraint = nvsite_restraint 00971 IF (PRESENT(g4x6_list)) molecule_kind%g4x6_list => g4x6_list 00972 IF (PRESENT(vsite_list)) molecule_kind%vsite_list => vsite_list 00973 IF (PRESENT(fixd_list)) molecule_kind%fixd_list => fixd_list 00974 IF (PRESENT(bond_kind_set)) molecule_kind%bond_kind_set => bond_kind_set 00975 IF (PRESENT(bend_kind_set)) molecule_kind%bend_kind_set => bend_kind_set 00976 IF (PRESENT(ub_kind_set)) molecule_kind%ub_kind_set => ub_kind_set 00977 IF (PRESENT(torsion_kind_set)) molecule_kind%torsion_kind_set => torsion_kind_set 00978 IF (PRESENT(impr_kind_set)) molecule_kind%impr_kind_set => impr_kind_set 00979 IF (PRESENT(opbend_kind_set)) molecule_kind%opbend_kind_set => opbend_kind_set 00980 IF (PRESENT(nshell)) molecule_kind%nshell = nshell 00981 IF (PRESENT(shell_list)) molecule_kind%shell_list => shell_list 00982 IF (PRESENT(molecule_list)) THEN 00983 n = SIZE(molecule_list) 00984 molecule_kind%nmolecule = n 00985 molecule_kind%molecule_list => molecule_list 00986 END IF 00987 00988 ELSE 00989 00990 CALL stop_program(routineN,moduleN,__LINE__,& 00991 "The pointer molecule_kind is not associated") 00992 00993 END IF 00994 00995 END SUBROUTINE set_molecule_kind 00996 00997 ! ***************************************************************************** 01003 SUBROUTINE write_molecule_kind(molecule_kind,output_unit,error) 01004 TYPE(molecule_kind_type), POINTER :: molecule_kind 01005 INTEGER, INTENT(in) :: output_unit 01006 TYPE(cp_error_type), INTENT(inout) :: error 01007 01008 CHARACTER(len=*), PARAMETER :: routineN = 'write_molecule_kind', 01009 routineP = moduleN//':'//routineN 01010 01011 CHARACTER(LEN=default_string_length) :: name 01012 INTEGER :: iatom, imolecule, natom, 01013 nmolecule 01014 LOGICAL :: failure 01015 TYPE(atomic_kind_type), POINTER :: atomic_kind 01016 01017 failure = .FALSE. 01018 IF (output_unit>0) THEN 01019 IF (ASSOCIATED(molecule_kind)) THEN 01020 natom = SIZE(molecule_kind%atom_list) 01021 nmolecule = SIZE(molecule_kind%molecule_list) 01022 01023 IF ( natom == 1 ) THEN 01024 atomic_kind => molecule_kind%atom_list(1)%atomic_kind 01025 CALL get_atomic_kind(atomic_kind=atomic_kind,name=name) 01026 WRITE (UNIT=output_unit,FMT="(/,T2,I5,A,T36,A,A,T64,A)")& 01027 molecule_kind%kind_number,& 01028 ". Molecule kind: "//TRIM(molecule_kind%name),& 01029 "Atomic kind name: ",TRIM(name) 01030 WRITE (UNIT=output_unit,FMT="(T9,A,L1,T55,A,T75,I6)")& 01031 "Automatic name: ",molecule_kind%molname_generated,& 01032 "Number of molecules:",nmolecule 01033 ELSE 01034 WRITE (UNIT=output_unit,FMT="(/,T2,I5,A,T50,A,T75,I6,/,T22,A)")& 01035 molecule_kind%kind_number,& 01036 ". Molecule kind: "//TRIM(molecule_kind%name),& 01037 "Number of atoms: ",natom,& 01038 "Atom Atomic kind name" 01039 DO iatom=1,natom 01040 atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind 01041 CALL get_atomic_kind(atomic_kind=atomic_kind,name=name) 01042 WRITE (UNIT=output_unit,FMT="(T20,I6,(7X,A18))")& 01043 iatom,TRIM(name) 01044 END DO 01045 WRITE (UNIT=output_unit,FMT="(/,T9,A,L1)")& 01046 "The name was automatically generated: ",& 01047 molecule_kind%molname_generated 01048 WRITE (UNIT=output_unit,FMT="(T9,A,I6,/,T9,A,(T30,5I10))")& 01049 "Number of molecules: ",nmolecule,"Molecule list:",& 01050 (molecule_kind%molecule_list(imolecule),imolecule=1,nmolecule) 01051 IF ( molecule_kind%nbond > 0 ) & 01052 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01053 "Number of bonds: ",molecule_kind%nbond 01054 IF ( molecule_kind%nbend > 0 ) & 01055 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01056 "Number of bends: ",molecule_kind%nbend 01057 IF ( molecule_kind%nub > 0 ) & 01058 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01059 "Number of Urey-Bradley:",molecule_kind%nub 01060 IF ( molecule_kind%ntorsion > 0 ) & 01061 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01062 "Number of torsions: ",molecule_kind%ntorsion 01063 IF ( molecule_kind%nimpr > 0 ) & 01064 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01065 "Number of improper: ",molecule_kind%nimpr 01066 IF ( molecule_kind%nopbend > 0 ) & 01067 WRITE (UNIT=output_unit,FMT="(1X,A30,I6)")& 01068 "Number of out opbends: ",molecule_kind%nopbend 01069 END IF 01070 01071 ELSE 01072 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01073 END IF 01074 01075 END IF 01076 END SUBROUTINE write_molecule_kind 01077 01078 ! ***************************************************************************** 01084 SUBROUTINE write_molecule_kind_set(molecule_kind_set,subsys_section,error) 01085 TYPE(molecule_kind_type), DIMENSION(:), 01086 POINTER :: molecule_kind_set 01087 TYPE(section_vals_type), POINTER :: subsys_section 01088 TYPE(cp_error_type), INTENT(inout) :: error 01089 01090 CHARACTER(len=*), PARAMETER :: routineN = 'write_molecule_kind_set', 01091 routineP = moduleN//':'//routineN 01092 01093 INTEGER :: handle, imolecule_kind, natom, nbend, nbond, nimpr, nmolecule, 01094 nmolecule_kind, nopbend, ntors, ntotal, nub, output_unit 01095 LOGICAL :: all_single_atoms, failure 01096 TYPE(cp_logger_type), POINTER :: logger 01097 TYPE(molecule_kind_type), POINTER :: molecule_kind 01098 01099 CALL timeset(routineN,handle) 01100 01101 failure = .FALSE. 01102 NULLIFY(logger) 01103 logger => cp_error_get_logger(error) 01104 output_unit = cp_print_key_unit_nr(logger,subsys_section,& 01105 "PRINT%MOLECULES",extension=".Log",error=error) 01106 IF (output_unit>0) THEN 01107 IF (ASSOCIATED(molecule_kind_set)) THEN 01108 WRITE (UNIT=output_unit,FMT="(/,/,T2,A)") "MOLECULE KIND INFORMATION" 01109 01110 nmolecule_kind = SIZE(molecule_kind_set) 01111 01112 all_single_atoms=.TRUE. 01113 DO imolecule_kind=1,nmolecule_kind 01114 natom = SIZE(molecule_kind_set(imolecule_kind)%atom_list) 01115 nmolecule = SIZE(molecule_kind_set(imolecule_kind)%molecule_list) 01116 IF (natom*nmolecule>1) all_single_atoms=.FALSE. 01117 ENDDO 01118 01119 IF (all_single_atoms) THEN 01120 WRITE (UNIT=output_unit,FMT="(/,/,T2,A)") & 01121 "All atoms are their own molecule, skipping detailed information" 01122 ELSE 01123 DO imolecule_kind=1,nmolecule_kind 01124 molecule_kind => molecule_kind_set(imolecule_kind) 01125 CALL write_molecule_kind(molecule_kind,output_unit,error) 01126 END DO 01127 ENDIF 01128 01129 CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set,& 01130 nbond=nbond,& 01131 nbend=nbend,& 01132 nub=nub,& 01133 ntorsion=ntors,& 01134 nimpr=nimpr,& 01135 nopbend=nopbend) 01136 ntotal = nbond+nbend+nub+ntors+nimpr+nopbend 01137 IF ( ntotal > 0 ) THEN 01138 WRITE(UNIT=output_unit,FMT="(/,/,T2,A,T45,A30,I6)") & 01139 "MOLECULE KIND SET INFORMATION", & 01140 "Total Number of bonds: ",nbond 01141 WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& 01142 "Total Number of bends: ",nbend 01143 WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& 01144 "Total Number of Urey-Bradley:",nub 01145 WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& 01146 "Total Number of torsions: ",ntors 01147 WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& 01148 "Total Number of improper: ",nimpr 01149 WRITE (UNIT=output_unit,FMT="(T45,A30,I6)")& 01150 "Total Number of opbends: ",nopbend 01151 END IF 01152 ELSE 01153 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01154 END IF 01155 01156 END IF 01157 CALL cp_print_key_finished_output(output_unit,logger,subsys_section,& 01158 "PRINT%MOLECULES",error=error) 01159 01160 CALL timestop(handle) 01161 01162 END SUBROUTINE write_molecule_kind_set 01163 01164 END MODULE molecule_kind_types
1.7.3