|
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 ! ***************************************************************************** 00017 MODULE force_fields_util 00018 USE atomic_kind_types, ONLY: atomic_kind_type,& 00019 get_atomic_kind 00020 USE cell_types, ONLY: cell_type 00021 USE colvar_types, ONLY: dist_colvar_id 00022 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00023 cp_print_key_unit_nr 00024 USE cp_units, ONLY: cp_unit_to_cp2k 00025 USE ewald_environment_types, ONLY: ewald_env_get,& 00026 ewald_environment_type 00027 USE f77_blas 00028 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_create,& 00029 fist_nonbond_env_set,& 00030 fist_nonbond_env_type 00031 USE force_field_types, ONLY: & 00032 allocate_bend_kind_set, allocate_bond_kind_set, & 00033 allocate_impr_kind_set, allocate_opbend_kind_set, & 00034 allocate_torsion_kind_set, allocate_ub_kind_set, amber_info_type, & 00035 bend_kind_type, bond_kind_type, charmm_info_type, & 00036 deallocate_bend_kind_set, deallocate_bond_kind_set, force_field_type, & 00037 gromos_info_type, impr_kind_dealloc_ref, impr_kind_type, & 00038 input_info_type, opbend_kind_type, torsion_kind_dealloc_ref, & 00039 torsion_kind_type, ub_kind_dealloc_ref, ub_kind_type 00040 USE force_fields_all, ONLY: & 00041 force_field_pack_bend, force_field_pack_bond, force_field_pack_charge, & 00042 force_field_pack_charges, force_field_pack_damp, & 00043 force_field_pack_eicut, force_field_pack_impr, & 00044 force_field_pack_nonbond, force_field_pack_nonbond14, & 00045 force_field_pack_opbend, force_field_pack_pol, & 00046 force_field_pack_radius, force_field_pack_shell, & 00047 force_field_pack_splines, force_field_pack_tors, force_field_pack_ub, & 00048 force_field_unique_bend, force_field_unique_bond, & 00049 force_field_unique_impr, force_field_unique_opbend, & 00050 force_field_unique_tors, force_field_unique_ub 00051 USE input_constants, ONLY: do_ff_undef 00052 USE input_section_types, ONLY: section_vals_get,& 00053 section_vals_get_subs_vals,& 00054 section_vals_type,& 00055 section_vals_val_get 00056 USE kinds, ONLY: default_path_length,& 00057 default_string_length,& 00058 dp 00059 USE memory_utilities, ONLY: reallocate 00060 USE molecule_kind_types, ONLY: & 00061 atom_type, bend_type, bond_type, colvar_constraint_type, & 00062 g3x3_constraint_type, g4x6_constraint_type, get_molecule_kind, & 00063 impr_type, molecule_kind_type, opbend_type, set_molecule_kind, & 00064 torsion_type, ub_type 00065 USE molecule_types_new, ONLY: get_molecule,& 00066 molecule_type 00067 USE pair_potential_types, ONLY: pair_potential_pp_type 00068 USE particle_types, ONLY: particle_type 00069 USE qmmm_types, ONLY: qmmm_env_mm_type 00070 USE shell_potential_types, ONLY: shell_kind_type 00071 USE string_utilities, ONLY: compress 00072 USE termination, ONLY: stop_program 00073 USE timings, ONLY: timeset,& 00074 timestop 00075 #include "cp_common_uses.h" 00076 00077 IMPLICIT NONE 00078 00079 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_util' 00080 00081 PRIVATE 00082 PUBLIC :: force_field_pack,& 00083 force_field_qeff_output,& 00084 clean_intra_force_kind,& 00085 get_generic_info 00086 00087 CONTAINS 00088 00089 ! ***************************************************************************** 00092 SUBROUTINE force_field_pack (particle_set,atomic_kind_set,molecule_kind_set, & 00093 molecule_set, ewald_env,fist_nonbond_env,ff_type,root_section, qmmm, & 00094 qmmm_env, mm_section, subsys_section, shell_particle_set,core_particle_set,& 00095 cell, error) 00096 00097 TYPE(particle_type), DIMENSION(:), 00098 POINTER :: particle_set 00099 TYPE(atomic_kind_type), DIMENSION(:), 00100 POINTER :: atomic_kind_set 00101 TYPE(molecule_kind_type), DIMENSION(:), 00102 POINTER :: molecule_kind_set 00103 TYPE(molecule_type), DIMENSION(:), 00104 POINTER :: molecule_set 00105 TYPE(ewald_environment_type), POINTER :: ewald_env 00106 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 00107 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00108 TYPE(section_vals_type), POINTER :: root_section 00109 LOGICAL, INTENT(IN), OPTIONAL :: qmmm 00110 TYPE(qmmm_env_mm_type), OPTIONAL, 00111 POINTER :: qmmm_env 00112 TYPE(section_vals_type), POINTER :: mm_section, subsys_section 00113 TYPE(particle_type), DIMENSION(:), 00114 POINTER :: shell_particle_set, 00115 core_particle_set 00116 TYPE(cell_type), POINTER :: cell 00117 TYPE(cp_error_type), INTENT(INOUT) :: error 00118 00119 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack', 00120 routineP = moduleN//':'//routineN 00121 00122 CHARACTER(LEN=default_string_length), 00123 DIMENSION(:), POINTER :: Ainfo 00124 INTEGER :: handle, iw, iw2, iw3, iw4, 00125 output_unit 00126 LOGICAL :: do_zbl, explicit, failure, 00127 fatal, ignore_fatal, my_qmmm 00128 REAL(KIND=dp) :: ewald_rcut, verlet_skin 00129 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 00130 TYPE(amber_info_type), POINTER :: amb_info 00131 TYPE(charmm_info_type), POINTER :: chm_info 00132 TYPE(cp_logger_type), POINTER :: logger 00133 TYPE(gromos_info_type), POINTER :: gro_info 00134 TYPE(input_info_type), POINTER :: inp_info 00135 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond, 00136 potparm_nonbond14 00137 TYPE(section_vals_type), POINTER :: charges_section 00138 00139 CALL timeset(routineN,handle) 00140 failure = .FALSE. 00141 fatal = .FALSE. 00142 ignore_fatal = ff_type%ignore_missing_critical 00143 NULLIFY(logger, Ainfo, charges_section, charges) 00144 logger => cp_error_get_logger(error) 00145 ! Error unit 00146 output_unit= cp_logger_get_default_io_unit(logger) 00147 00148 iw = cp_print_key_unit_nr(logger,mm_section,"PRINT%FF_INFO",& 00149 extension=".mmLog",error=error) 00150 iw2= cp_print_key_unit_nr(logger,mm_section,"PRINT%FF_INFO/SPLINE_INFO",& 00151 extension=".mmLog",error=error) 00152 iw3= cp_print_key_unit_nr(logger,mm_section,"PRINT%FF_INFO/SPLINE_DATA",& 00153 extension=".mmLog",error=error) 00154 iw4= cp_print_key_unit_nr(logger,mm_section,"PRINT%PROGRAM_RUN_INFO",& 00155 extension=".mmLog",error=error) 00156 NULLIFY(potparm_nonbond14, potparm_nonbond) 00157 my_qmmm = .FALSE. 00158 IF (PRESENT(qmmm).AND.PRESENT(qmmm_env)) my_qmmm = qmmm 00159 inp_info => ff_type%inp_info 00160 chm_info => ff_type%chm_info 00161 gro_info => ff_type%gro_info 00162 amb_info => ff_type%amb_info 00163 !----------------------------------------------------------------------------- 00164 !----------------------------------------------------------------------------- 00165 ! 1. Determine the number of unique bond kind and allocate bond_kind_set 00166 !----------------------------------------------------------------------------- 00167 CALL force_field_unique_bond (particle_set, molecule_kind_set, molecule_set, & 00168 ff_type, error) 00169 !----------------------------------------------------------------------------- 00170 !----------------------------------------------------------------------------- 00171 ! 2. Determine the number of unique bend kind and allocate bend_kind_set 00172 !----------------------------------------------------------------------------- 00173 CALL force_field_unique_bend (particle_set, molecule_kind_set, molecule_set, & 00174 ff_type, error) 00175 !----------------------------------------------------------------------------- 00176 !----------------------------------------------------------------------------- 00177 ! 3. Determine the number of unique Urey-Bradley kind and allocate ub_kind_set 00178 !----------------------------------------------------------------------------- 00179 CALL force_field_unique_ub (particle_set, molecule_kind_set, molecule_set, & 00180 error) 00181 !----------------------------------------------------------------------------- 00182 !----------------------------------------------------------------------------- 00183 ! 4. Determine the number of unique torsion kind and allocate torsion_kind_set 00184 !----------------------------------------------------------------------------- 00185 CALL force_field_unique_tors (particle_set, molecule_kind_set, molecule_set, & 00186 ff_type, error) 00187 !----------------------------------------------------------------------------- 00188 !----------------------------------------------------------------------------- 00189 ! 5. Determine the number of unique impr kind and allocate impr_kind_set 00190 !----------------------------------------------------------------------------- 00191 CALL force_field_unique_impr (particle_set, molecule_kind_set, molecule_set, & 00192 ff_type, error) 00193 !----------------------------------------------------------------------------- 00194 !----------------------------------------------------------------------------- 00195 ! 6. Determine the number of unique opbend kind and allocate opbend_kind_set 00196 !----------------------------------------------------------------------------- 00197 CALL force_field_unique_opbend (particle_set, molecule_kind_set, molecule_set, & 00198 ff_type, error) 00199 !----------------------------------------------------------------------------- 00200 !----------------------------------------------------------------------------- 00201 ! 7. Bonds 00202 !----------------------------------------------------------------------------- 00203 CALL force_field_pack_bond (particle_set, molecule_kind_set, molecule_set, & 00204 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) 00205 !----------------------------------------------------------------------------- 00206 !----------------------------------------------------------------------------- 00207 ! 8. Bends 00208 !----------------------------------------------------------------------------- 00209 CALL force_field_pack_bend (particle_set, molecule_kind_set, molecule_set, & 00210 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) 00211 ! Give information and abort if any bond or angle parameter is missing.. 00212 CALL release_FF_missing_par(fatal,ignore_fatal,AInfo,output_unit,iw,error) 00213 !----------------------------------------------------------------------------- 00214 !----------------------------------------------------------------------------- 00215 ! 9. Urey-Bradley 00216 !----------------------------------------------------------------------------- 00217 CALL force_field_pack_ub (particle_set, molecule_kind_set, molecule_set, & 00218 Ainfo, chm_info, inp_info, iw, error) 00219 !----------------------------------------------------------------------------- 00220 !----------------------------------------------------------------------------- 00221 ! 10. Torsions 00222 !----------------------------------------------------------------------------- 00223 CALL force_field_pack_tors (particle_set, molecule_kind_set, molecule_set, & 00224 Ainfo, chm_info, inp_info, gro_info, amb_info, iw, error) 00225 !----------------------------------------------------------------------------- 00226 !----------------------------------------------------------------------------- 00227 ! 11. Impropers 00228 !----------------------------------------------------------------------------- 00229 CALL force_field_pack_impr (particle_set, molecule_kind_set, molecule_set, & 00230 Ainfo, chm_info, inp_info, gro_info, amb_info, error) 00231 !----------------------------------------------------------------------------- 00232 !----------------------------------------------------------------------------- 00233 ! 12. Out of plane bends 00234 !----------------------------------------------------------------------------- 00235 CALL force_field_pack_opbend (particle_set, molecule_kind_set, molecule_set, & 00236 Ainfo, inp_info, error) 00237 ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing 00238 ! continue calculation.. 00239 CALL release_FF_missing_par(fatal,ignore_fatal,AInfo,output_unit,iw,error) 00240 00241 charges_section => section_vals_get_subs_vals(mm_section,"FORCEFIELD%CHARGES",error=error) 00242 CALL section_vals_get(charges_section, explicit=explicit, error=error) 00243 IF (.NOT.explicit) THEN 00244 !----------------------------------------------------------------------------- 00245 !----------------------------------------------------------------------------- 00246 ! 13.a Set up atomic_kind_set()%fist_potentail%[qeff] and shell 00247 ! potential parameters 00248 !----------------------------------------------------------------------------- 00249 CALL force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4,& 00250 Ainfo, my_qmmm, inp_info, error) 00251 ! Give information only if charge is missing and abort.. 00252 CALL release_FF_missing_par(fatal,ignore_fatal,AInfo,output_unit,iw,error) 00253 ELSE 00254 !----------------------------------------------------------------------------- 00255 !----------------------------------------------------------------------------- 00256 ! 13.b Setup a global array of classical charges - this avoids the packing and 00257 ! allows the usage of different charges for same atomic types 00258 !----------------------------------------------------------------------------- 00259 CALL force_field_pack_charges(charges, charges_section, particle_set, my_qmmm,& 00260 qmmm_env, inp_info, iw4, error) 00261 END IF 00262 00263 !----------------------------------------------------------------------------- 00264 !----------------------------------------------------------------------------- 00265 ! 14. Set up the radius of the electrostatic multipole in Fist 00266 !----------------------------------------------------------------------------- 00267 CALL force_field_pack_radius(atomic_kind_set, iw, subsys_section, error) 00268 00269 !----------------------------------------------------------------------------- 00270 !----------------------------------------------------------------------------- 00271 ! 15. Set up the polarizable FF parameters 00272 !----------------------------------------------------------------------------- 00273 CALL force_field_pack_pol(atomic_kind_set, iw, inp_info, error) 00274 CALL force_field_pack_damp(atomic_kind_set, iw, inp_info, error) 00275 00276 !----------------------------------------------------------------------------- 00277 !----------------------------------------------------------------------------- 00278 ! 16. Set up Shell potential 00279 !----------------------------------------------------------------------------- 00280 CALL force_field_pack_shell ( particle_set, atomic_kind_set,& 00281 molecule_kind_set, molecule_set, root_section, subsys_section,& 00282 shell_particle_set, core_particle_set, cell, iw, inp_info, error) 00283 IF (ff_type%do_nonbonded) THEN 00284 !----------------------------------------------------------------------------- 00285 !----------------------------------------------------------------------------- 00286 ! 17. Set up potparm_nonbond14 00287 !----------------------------------------------------------------------------- 00288 ! Move the data from the info structures to potparm_nonbond 00289 CALL force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, Ainfo,& 00290 chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env, & 00291 iw2, iw3, error) 00292 ! Give information if any 1-4 is missing.. continue calculation.. 00293 CALL release_FF_missing_par(fatal,ignore_fatal,AInfo,output_unit,iw,error) 00294 ! Create the spline data 00295 CALL section_vals_val_get(mm_section,"FORCEFIELD%ZBL_SCATTERING",l_val=do_zbl,error=error) 00296 CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, & 00297 potparm_nonbond14, do_zbl, nonbonded_type="NONBONDED14", error=error) 00298 !----------------------------------------------------------------------------- 00299 !----------------------------------------------------------------------------- 00300 ! 18. Set up potparm_nonbond 00301 !----------------------------------------------------------------------------- 00302 ! Move the data from the info structures to potparm_nonbond 00303 CALL force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, & 00304 fatal, iw, Ainfo, chm_info, inp_info, gro_info, amb_info, & 00305 potparm_nonbond, ewald_env, error) 00306 ! Give information and abort if any pair potential spline is missing.. 00307 CALL release_FF_missing_par(fatal,ignore_fatal,AInfo,output_unit,iw,error) 00308 ! Create the spline data 00309 CALL section_vals_val_get(mm_section,"FORCEFIELD%ZBL_SCATTERING",l_val=do_zbl,error=error) 00310 CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, & 00311 potparm_nonbond, do_zbl, nonbonded_type="NONBONDED", error=error) 00312 END IF 00313 !----------------------------------------------------------------------------- 00314 !----------------------------------------------------------------------------- 00315 ! 19. Create nonbond environment 00316 !----------------------------------------------------------------------------- 00317 CALL ewald_env_get(ewald_env, rcut=ewald_rcut, error=error) 00318 CALL section_vals_val_get(mm_section,"NEIGHBOR_LISTS%VERLET_SKIN",& 00319 r_val=verlet_skin,error=error) 00320 CALL fist_nonbond_env_create (fist_nonbond_env, atomic_kind_set, & 00321 potparm_nonbond14, potparm_nonbond, ff_type%do_nonbonded, & 00322 verlet_skin, ewald_rcut, ff_type%ei_scale14, & 00323 ff_type%vdw_scale14, ff_type%shift_cutoff, error) 00324 CALL fist_nonbond_env_set(fist_nonbond_env, charges=charges, error=error) 00325 ! Compute the electrostatic interaction cutoffs. 00326 CALL force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, & 00327 ewald_env, error) 00328 00329 CALL cp_print_key_finished_output(iw4,logger,mm_section,"PRINT%PROGRAM_RUN_INFO",error=error) 00330 CALL cp_print_key_finished_output(iw3,logger,mm_section,"PRINT%FF_INFO/SPLINE_DATA",error=error) 00331 CALL cp_print_key_finished_output(iw2,logger,mm_section,"PRINT%FF_INFO/SPLINE_INFO",error=error) 00332 CALL cp_print_key_finished_output(iw,logger,mm_section,"PRINT%FF_INFO",error=error) 00333 CALL timestop(handle) 00334 00335 END SUBROUTINE force_field_pack 00336 00337 ! ***************************************************************************** 00340 SUBROUTINE release_FF_missing_par(fatal, ignore_fatal, array, output_unit, iw, error) 00341 LOGICAL, INTENT(INOUT) :: fatal, ignore_fatal 00342 CHARACTER(LEN=default_string_length), 00343 DIMENSION(:), POINTER :: array 00344 INTEGER, INTENT(IN) :: output_unit, iw 00345 TYPE(cp_error_type), INTENT(inout) :: error 00346 00347 CHARACTER(len=*), PARAMETER :: routineN = 'release_FF_missing_par', 00348 routineP = moduleN//':'//routineN 00349 00350 INTEGER :: i, stat 00351 LOGICAL :: failure 00352 00353 failure = .FALSE. 00354 IF (ASSOCIATED(array)) THEN 00355 IF (output_unit>0) THEN 00356 WRITE(output_unit,*) 00357 WRITE(output_unit,'(T2,"FORCEFIELD|",A)')& 00358 " WARNING: A non Critical ForceField parameter is missing! CP2K GOES!",& 00359 " Non critical parameters are:Urey-Bradley,Torsions,Impropers, Opbends and 1-4",& 00360 " All missing parameters will not contribute to the potential energy!" 00361 IF (fatal.OR.iw>0) THEN 00362 WRITE(output_unit,*) 00363 DO i = 1,SIZE(array) 00364 WRITE(output_unit,'(A)')array(i) 00365 END DO 00366 END IF 00367 IF (.NOT.fatal.AND.iw<0) THEN 00368 WRITE(output_unit,'(T2,"FORCEFIELD|",A)')& 00369 " Activate the print key FF_INFO to have a list of missing parameters" 00370 WRITE(output_unit,*) 00371 END IF 00372 END IF 00373 DEALLOCATE(array,stat=stat) 00374 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00375 END IF 00376 IF (fatal) THEN 00377 IF (ignore_fatal) THEN 00378 IF (output_unit>0) THEN 00379 WRITE(output_unit,*) 00380 WRITE(output_unit,'(T2,"FORCEFIELD|",A)')& 00381 " WARNING: Ignoring missing critical FF parameters! CP2K GOES!",& 00382 " Critical parameters are: Bonds, Bends and Charges",& 00383 " All missing parameters will not contribute to the potential energy!" 00384 END IF 00385 ELSE 00386 CALL stop_program(routineN,moduleN,__LINE__,& 00387 "Missing critical ForceField parameters!") 00388 END IF 00389 END IF 00390 END SUBROUTINE release_FF_missing_par 00391 00392 ! ***************************************************************************** 00395 SUBROUTINE force_field_qeff_output (particle_set,molecule_kind_set,& 00396 molecule_set,mm_section,charges,error) 00397 00398 TYPE(particle_type), DIMENSION(:), 00399 POINTER :: particle_set 00400 TYPE(molecule_kind_type), DIMENSION(:), 00401 POINTER :: molecule_kind_set 00402 TYPE(molecule_type), DIMENSION(:), 00403 POINTER :: molecule_set 00404 TYPE(section_vals_type), POINTER :: mm_section 00405 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 00406 TYPE(cp_error_type), INTENT(inout) :: error 00407 00408 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_qeff_output', 00409 routineP = moduleN//':'//routineN 00410 00411 CHARACTER(LEN=default_string_length) :: atmname, molname 00412 INTEGER :: first, handle, iatom, imol, 00413 iw, j, jatom 00414 LOGICAL :: failure, shell_active 00415 REAL(KIND=dp) :: mass, mass_mol, mass_sum, 00416 qeff, qeff_mol, qeff_sum 00417 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 00418 TYPE(atomic_kind_type), POINTER :: atomic_kind 00419 TYPE(cp_logger_type), POINTER :: logger 00420 TYPE(molecule_kind_type), POINTER :: molecule_kind 00421 TYPE(molecule_type), POINTER :: molecule 00422 TYPE(shell_kind_type), POINTER :: shell 00423 00424 CALL timeset(routineN,handle) 00425 failure = .FALSE. 00426 NULLIFY(logger) 00427 logger => cp_error_get_logger(error) 00428 iw = cp_print_key_unit_nr(logger,mm_section,"PRINT%FF_INFO",& 00429 extension=".mmLog",error=error) 00430 00431 qeff = 0.0_dp 00432 qeff_mol = 0.0_dp 00433 qeff_sum = 0.0_dp 00434 mass_sum = 0.0_dp 00435 00436 !----------------------------------------------------------------------------- 00437 !----------------------------------------------------------------------------- 00438 ! 1. Sum of [qeff,mass] for each molecule_kind 00439 !----------------------------------------------------------------------------- 00440 DO imol=1,SIZE(molecule_kind_set) 00441 qeff_mol=0.0_dp 00442 mass_mol=0.0_dp 00443 molecule_kind => molecule_kind_set(imol) 00444 00445 j=molecule_kind_set(imol)%molecule_list(1) 00446 molecule => molecule_set(j) 00447 CALL get_molecule(molecule=molecule, first_atom=first) 00448 00449 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00450 name=molname,atom_list=atom_list) 00451 DO iatom=1,SIZE(atom_list) 00452 atomic_kind => atom_list(iatom)%atomic_kind 00453 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00454 name=atmname,qeff=qeff,mass=mass, shell_active=shell_active, shell=shell) 00455 IF(shell_active) qeff = shell%charge_core+shell%charge_shell 00456 IF (ASSOCIATED(charges)) THEN 00457 jatom=first-1+iatom 00458 qeff=charges(jatom) 00459 END IF 00460 IF(iw>0) WRITE(iw,*) " atom ",iatom," ",TRIM(atmname)," charge = ",qeff 00461 qeff_mol = qeff_mol + qeff 00462 mass_mol = mass_mol + mass 00463 END DO 00464 CALL set_molecule_kind(molecule_kind=molecule_kind,charge=qeff_mol,mass=mass_mol) 00465 IF(iw>0) WRITE(iw,*) " Mol Kind ",TRIM(molname)," charge = ",qeff_mol 00466 END DO 00467 00468 !----------------------------------------------------------------------------- 00469 !----------------------------------------------------------------------------- 00470 ! 2. Sum of [qeff,mass] for particle_set 00471 !----------------------------------------------------------------------------- 00472 DO iatom=1,SIZE(particle_set) 00473 atomic_kind => particle_set(iatom)%atomic_kind 00474 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00475 name=atmname,qeff=qeff,mass=mass,shell_active=shell_active, shell=shell) 00476 IF(shell_active) qeff = shell%charge_core+shell%charge_shell 00477 IF (ASSOCIATED(charges)) THEN 00478 qeff=charges(iatom) 00479 END IF 00480 IF(iw>0) WRITE(iw,*) " atom ",iatom," ",TRIM(atmname),& 00481 " charge = ",qeff 00482 qeff_sum = qeff_sum + qeff 00483 mass_sum = mass_sum + mass 00484 END DO 00485 IF(iw>0) WRITE(iw,'(A,F20.10)') " Total system charge = ",qeff_sum 00486 IF(iw>0) WRITE(iw,'(A,F20.10)') " Total system mass = ",mass_sum 00487 00488 CALL cp_print_key_finished_output(iw,logger,mm_section,& 00489 "PRINT%FF_INFO",error=error) 00490 CALL timestop(handle) 00491 END SUBROUTINE force_field_qeff_output 00492 00493 ! ***************************************************************************** 00496 SUBROUTINE clean_intra_force_kind (molecule_kind_set,mm_section,error) 00497 00498 TYPE(molecule_kind_type), DIMENSION(:), 00499 POINTER :: molecule_kind_set 00500 TYPE(section_vals_type), POINTER :: mm_section 00501 TYPE(cp_error_type), INTENT(inout) :: error 00502 00503 CHARACTER(len=*), PARAMETER :: routineN = 'clean_intra_force_kind', 00504 routineP = moduleN//':'//routineN 00505 00506 INTEGER :: atm2_a, atm2_b, atm2_c, atm_a, atm_b, atm_c, atm_d, counter, 00507 handle, i, ibend, ibond, icolv, ig3x3, ig4x6, iimpr, ikind, iopbend, 00508 itorsion, iub, iw, j, k, nbend, nbond, newkind, ng3x3, ng4x6, nimpr, 00509 nopbend, ntorsion, nub, stat 00510 INTEGER, POINTER :: bad1(:), bad2(:) 00511 LOGICAL :: failure, unsetme, valid_kind 00512 TYPE(bend_kind_type), DIMENSION(:), 00513 POINTER :: bend_kind_set, 00514 new_bend_kind_set 00515 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list, new_bend_list 00516 TYPE(bond_kind_type), DIMENSION(:), 00517 POINTER :: bond_kind_set, 00518 new_bond_kind_set 00519 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list, new_bond_list 00520 TYPE(colvar_constraint_type), 00521 DIMENSION(:), POINTER :: colv_list 00522 TYPE(cp_logger_type), POINTER :: logger 00523 TYPE(g3x3_constraint_type), 00524 DIMENSION(:), POINTER :: g3x3_list 00525 TYPE(g4x6_constraint_type), 00526 DIMENSION(:), POINTER :: g4x6_list 00527 TYPE(impr_kind_type), DIMENSION(:), 00528 POINTER :: impr_kind_set, 00529 new_impr_kind_set 00530 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list, new_impr_list 00531 TYPE(molecule_kind_type), POINTER :: molecule_kind 00532 TYPE(opbend_kind_type), DIMENSION(:), 00533 POINTER :: new_opbend_kind_set, 00534 opbend_kind_set 00535 TYPE(opbend_type), DIMENSION(:), POINTER :: new_opbend_list, opbend_list 00536 TYPE(torsion_kind_type), DIMENSION(:), 00537 POINTER :: new_torsion_kind_set, 00538 torsion_kind_set 00539 TYPE(torsion_type), DIMENSION(:), 00540 POINTER :: new_torsion_list, torsion_list 00541 TYPE(ub_kind_type), DIMENSION(:), 00542 POINTER :: new_ub_kind_set, ub_kind_set 00543 TYPE(ub_type), DIMENSION(:), POINTER :: new_ub_list, ub_list 00544 00545 CALL timeset(routineN,handle) 00546 failure = .FALSE. 00547 NULLIFY(logger) 00548 logger => cp_error_get_logger(error) 00549 iw = cp_print_key_unit_nr(logger,mm_section,"PRINT%FF_INFO",& 00550 extension=".mmLog",error=error) 00551 !----------------------------------------------------------------------------- 00552 !----------------------------------------------------------------------------- 00553 ! 1. Lets Tag the unwanted bonds due to the use of distance constraint 00554 !----------------------------------------------------------------------------- 00555 DO ikind=1,SIZE(molecule_kind_set) 00556 molecule_kind => molecule_kind_set(ikind) 00557 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00558 colv_list=colv_list,& 00559 nbond=nbond,& 00560 bond_list=bond_list) 00561 IF (ASSOCIATED(colv_list)) THEN 00562 DO icolv=1,SIZE(colv_list) 00563 IF ((colv_list(icolv)%type_id == dist_colvar_id).AND.& 00564 ((.NOT.colv_list(icolv)%use_points).OR.(SIZE(colv_list(icolv)%i_atoms)==2))) THEN 00565 atm_a = colv_list(icolv)%i_atoms(1) 00566 atm_b = colv_list(icolv)%i_atoms(2) 00567 DO ibond=1,nbond 00568 unsetme = .FALSE. 00569 atm2_a = bond_list(ibond)%a 00570 atm2_b = bond_list(ibond)%b 00571 IF(atm2_a==atm_a .AND. atm2_b==atm_b) unsetme=.TRUE. 00572 IF(atm2_a==atm_b .AND. atm2_b==atm_a) unsetme=.TRUE. 00573 IF(unsetme) bond_list(ibond)%id_type = do_ff_undef 00574 END DO 00575 END IF 00576 END DO 00577 END IF 00578 END DO 00579 !----------------------------------------------------------------------------- 00580 !----------------------------------------------------------------------------- 00581 ! 2. Lets Tag the unwanted bends due to the use of distance constraint 00582 !----------------------------------------------------------------------------- 00583 DO ikind=1,SIZE(molecule_kind_set) 00584 molecule_kind => molecule_kind_set(ikind) 00585 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00586 colv_list=colv_list,& 00587 nbend=nbend,& 00588 bend_list=bend_list) 00589 IF (ASSOCIATED(colv_list)) THEN 00590 DO ibend=1,nbend 00591 unsetme = .FALSE. 00592 atm_a = bend_list(ibend)%a 00593 atm_b = bend_list(ibend)%b 00594 atm_c = bend_list(ibend)%c 00595 DO icolv=1,SIZE(colv_list) 00596 IF ((colv_list(icolv)%type_id == dist_colvar_id).AND.& 00597 ((.NOT.colv_list(icolv)%use_points).OR.(SIZE(colv_list(icolv)%i_atoms)==2))) THEN 00598 atm2_a = colv_list(icolv)%i_atoms(1) 00599 atm2_b = colv_list(icolv)%i_atoms(2) 00600 ! Check that the bonds we're constraining does not involve atoms defining 00601 ! a bend.. 00602 IF ((( atm2_a == atm_a ).AND.( atm2_b == atm_c )).OR.& 00603 (( atm2_a == atm_c ).AND.( atm2_b == atm_a ))) THEN 00604 unsetme=.TRUE. 00605 EXIT 00606 END IF 00607 END IF 00608 END DO 00609 IF(unsetme) bend_list(ibend)%id_type = do_ff_undef 00610 END DO 00611 END IF 00612 END DO 00613 !----------------------------------------------------------------------------- 00614 !----------------------------------------------------------------------------- 00615 ! 3. Lets Tag the unwanted bonds due to the use of 3x3 00616 !----------------------------------------------------------------------------- 00617 DO ikind=1,SIZE(molecule_kind_set) 00618 molecule_kind => molecule_kind_set(ikind) 00619 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00620 ng3x3=ng3x3,& 00621 g3x3_list=g3x3_list,& 00622 nbond=nbond,& 00623 bond_list=bond_list) 00624 DO ig3x3=1,ng3x3 00625 atm_a = g3x3_list(ig3x3)%a 00626 atm_b = g3x3_list(ig3x3)%b 00627 atm_c = g3x3_list(ig3x3)%c 00628 DO ibond=1,nbond 00629 unsetme = .FALSE. 00630 atm2_a = bond_list(ibond)%a 00631 atm2_b = bond_list(ibond)%b 00632 IF(atm2_a==atm_a .AND. atm2_b==atm_b) unsetme=.TRUE. 00633 IF(atm2_a==atm_a .AND. atm2_b==atm_c) unsetme=.TRUE. 00634 IF(atm2_a==atm_c .AND. atm2_b==atm_c) unsetme=.TRUE. 00635 IF(unsetme) bond_list(ibond)%id_type = do_ff_undef 00636 END DO 00637 END DO 00638 END DO 00639 !----------------------------------------------------------------------------- 00640 !----------------------------------------------------------------------------- 00641 ! 4. Lets Tag the unwanted bends due to the use of 3x3 00642 !----------------------------------------------------------------------------- 00643 DO ikind=1,SIZE(molecule_kind_set) 00644 molecule_kind => molecule_kind_set(ikind) 00645 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00646 ng3x3=ng3x3,& 00647 g3x3_list=g3x3_list,& 00648 nbend=nbend,& 00649 bend_list=bend_list) 00650 DO ig3x3=1,ng3x3 00651 atm_a = g3x3_list(ig3x3)%a 00652 atm_b = g3x3_list(ig3x3)%b 00653 atm_c = g3x3_list(ig3x3)%c 00654 DO ibend=1,nbend 00655 unsetme = .FALSE. 00656 atm2_a = bend_list(ibend)%a 00657 atm2_b = bend_list(ibend)%b 00658 atm2_c = bend_list(ibend)%c 00659 IF(atm2_a==atm_a .AND. atm2_b==atm_b .AND. atm2_c==atm_c) unsetme=.TRUE. 00660 IF(atm2_a==atm_a .AND. atm2_b==atm_c .AND. atm2_c==atm_b) unsetme=.TRUE. 00661 IF(atm2_a==atm_b .AND. atm2_b==atm_a .AND. atm2_c==atm_c) unsetme=.TRUE. 00662 IF(atm2_a==atm_b .AND. atm2_b==atm_c .AND. atm2_c==atm_a) unsetme=.TRUE. 00663 IF(unsetme) bend_list(ibend)%id_type = do_ff_undef 00664 END DO 00665 END DO 00666 END DO 00667 !----------------------------------------------------------------------------- 00668 !----------------------------------------------------------------------------- 00669 ! 5. Lets Tag the unwanted bonds due to the use of 4x6 00670 !----------------------------------------------------------------------------- 00671 DO ikind=1,SIZE(molecule_kind_set) 00672 molecule_kind => molecule_kind_set(ikind) 00673 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00674 ng4x6=ng4x6,& 00675 g4x6_list=g4x6_list,& 00676 nbond=nbond,& 00677 bond_list=bond_list) 00678 DO ig4x6=1,ng4x6 00679 atm_a = g4x6_list(ig4x6)%a 00680 atm_b = g4x6_list(ig4x6)%b 00681 atm_c = g4x6_list(ig4x6)%c 00682 atm_d = g4x6_list(ig4x6)%d 00683 DO ibond=1,nbond 00684 unsetme = .FALSE. 00685 atm2_a = bond_list(ibond)%a 00686 atm2_b = bond_list(ibond)%b 00687 IF(atm2_a==atm_a .AND. atm2_b==atm_b) unsetme=.TRUE. 00688 IF(atm2_a==atm_a .AND. atm2_b==atm_c) unsetme=.TRUE. 00689 IF(atm2_a==atm_a .AND. atm2_b==atm_d) unsetme=.TRUE. 00690 IF(unsetme) bond_list(ibond)%id_type = do_ff_undef 00691 END DO 00692 END DO 00693 END DO 00694 !----------------------------------------------------------------------------- 00695 !----------------------------------------------------------------------------- 00696 ! 6. Lets Tag the unwanted bends due to the use of 4x6 00697 !----------------------------------------------------------------------------- 00698 DO ikind=1,SIZE(molecule_kind_set) 00699 molecule_kind => molecule_kind_set(ikind) 00700 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00701 ng4x6=ng4x6,& 00702 g4x6_list=g4x6_list,& 00703 nbend=nbend,& 00704 bend_list=bend_list) 00705 DO ig4x6=1,ng4x6 00706 atm_a = g4x6_list(ig4x6)%a 00707 atm_b = g4x6_list(ig4x6)%b 00708 atm_c = g4x6_list(ig4x6)%c 00709 atm_d = g4x6_list(ig4x6)%d 00710 DO ibend=1,nbend 00711 unsetme = .FALSE. 00712 atm2_a = bend_list(ibend)%a 00713 atm2_b = bend_list(ibend)%b 00714 atm2_c = bend_list(ibend)%c 00715 IF(atm2_a==atm_b .AND. atm2_b==atm_a .AND. atm2_c==atm_c) unsetme=.TRUE. 00716 IF(atm2_a==atm_b .AND. atm2_b==atm_a .AND. atm2_c==atm_d) unsetme=.TRUE. 00717 IF(atm2_a==atm_c .AND. atm2_b==atm_a .AND. atm2_c==atm_d) unsetme=.TRUE. 00718 IF(unsetme) bend_list(ibend)%id_type = do_ff_undef 00719 END DO 00720 END DO 00721 END DO 00722 !----------------------------------------------------------------------------- 00723 !----------------------------------------------------------------------------- 00724 ! 7. Count the number of UNSET bond kinds there are 00725 !----------------------------------------------------------------------------- 00726 DO ikind=1,SIZE(molecule_kind_set) 00727 molecule_kind => molecule_kind_set(ikind) 00728 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00729 nbond=nbond,& 00730 bond_kind_set=bond_kind_set,& 00731 bond_list=bond_list) 00732 IF(nbond>0) THEN 00733 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old BOND Count: ",& 00734 SIZE(bond_list),SIZE(bond_kind_set) 00735 IF (iw>0) WRITE(iw,'(2I6)')(bond_list(ibond)%a,bond_list(ibond)%b,ibond=1,SIZE(bond_list)) 00736 NULLIFY(bad1,bad2) 00737 ALLOCATE(bad1(SIZE(bond_kind_set)),STAT=stat) 00738 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00739 bad1(:)=0 00740 DO ibond=1,SIZE(bond_kind_set) 00741 unsetme=.FALSE. 00742 IF(bond_kind_set(ibond)%id_type==do_ff_undef) unsetme = .TRUE. 00743 valid_kind=.FALSE. 00744 DO i=1,SIZE(bond_list) 00745 IF(bond_list(i)%id_type/=do_ff_undef.AND.& 00746 bond_list(i)%bond_kind%kind_number==ibond) THEN 00747 valid_kind=.TRUE. 00748 EXIT 00749 END IF 00750 END DO 00751 IF(.NOT.valid_kind) unsetme = .TRUE. 00752 IF(unsetme) bad1(ibond) = 1 00753 END DO 00754 IF(SUM(bad1)/=0) THEN 00755 counter = SIZE(bond_kind_set)-SUM(bad1) 00756 CALL allocate_bond_kind_set(new_bond_kind_set,counter,error) 00757 counter=0 00758 DO ibond=1,SIZE(bond_kind_set) 00759 IF(bad1(ibond)==0) THEN 00760 counter=counter+1 00761 new_bond_kind_set(counter)= bond_kind_set(ibond) 00762 END IF 00763 END DO 00764 counter=0 00765 ALLOCATE(bad2(SIZE(bond_list)),STAT=stat) 00766 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00767 bad2(:)=0 00768 DO ibond=1,SIZE(bond_list) 00769 unsetme = .FALSE. 00770 IF(bond_list(ibond)%bond_kind%id_type==do_ff_undef) unsetme = .TRUE. 00771 IF(bond_list(ibond)%id_type==do_ff_undef) unsetme = .TRUE. 00772 IF(unsetme) bad2(ibond) = 1 00773 END DO 00774 IF(SUM(bad2)/=0) THEN 00775 counter = SIZE(bond_list)-SUM(bad2) 00776 ALLOCATE(new_bond_list(counter),STAT=stat) 00777 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00778 counter=0 00779 DO ibond=1,SIZE(bond_list) 00780 IF(bad2(ibond)==0) THEN 00781 counter=counter+1 00782 new_bond_list(counter) = bond_list(ibond) 00783 newkind = bond_list(ibond)%bond_kind%kind_number 00784 newkind = newkind - SUM(bad1(1:newkind)) 00785 new_bond_list(counter)%bond_kind => new_bond_kind_set(newkind) 00786 END IF 00787 END DO 00788 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00789 nbond=SIZE(new_bond_list),& 00790 bond_kind_set=new_bond_kind_set,& 00791 bond_list=new_bond_list) 00792 DO ibond=1,SIZE(new_bond_kind_set) 00793 new_bond_kind_set(ibond)%kind_number=ibond 00794 END DO 00795 END IF 00796 DEALLOCATE(bad2,STAT=stat) 00797 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00798 CALL deallocate_bond_kind_set(bond_kind_set,error) 00799 DEALLOCATE(bond_list,STAT=stat) 00800 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00801 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New BOND Count: ",& 00802 SIZE(new_bond_list),SIZE(new_bond_kind_set) 00803 IF (iw>0) WRITE(iw,'(2I6)')(new_bond_list(ibond)%a,new_bond_list(ibond)%b,& 00804 ibond=1,SIZE(new_bond_list)) 00805 END IF 00806 DEALLOCATE(bad1,STAT=stat) 00807 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00808 END IF 00809 END DO 00810 !----------------------------------------------------------------------------- 00811 !----------------------------------------------------------------------------- 00812 ! 8. Count the number of UNSET bend kinds there are 00813 !----------------------------------------------------------------------------- 00814 DO ikind=1,SIZE(molecule_kind_set) 00815 molecule_kind => molecule_kind_set(ikind) 00816 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00817 nbend=nbend,& 00818 bend_kind_set=bend_kind_set,& 00819 bend_list=bend_list) 00820 IF(nbend>0) THEN 00821 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old BEND Count: ",& 00822 SIZE(bend_list),SIZE(bend_kind_set) 00823 IF (iw>0) WRITE(iw,'(3I6)')(bend_list(ibend)%a,bend_list(ibend)%b,& 00824 bend_list(ibend)%c,ibend=1,SIZE(bend_list)) 00825 NULLIFY(bad1,bad2) 00826 ALLOCATE(bad1(SIZE(bend_kind_set)),STAT=stat) 00827 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00828 bad1(:)=0 00829 DO ibend=1,SIZE(bend_kind_set) 00830 unsetme=.FALSE. 00831 IF(bend_kind_set(ibend)%id_type==do_ff_undef) unsetme = .TRUE. 00832 valid_kind=.FALSE. 00833 DO i=1,SIZE(bend_list) 00834 IF(bend_list(i)%id_type/=do_ff_undef.AND.& 00835 bend_list(i)%bend_kind%kind_number==ibend) THEN 00836 valid_kind=.TRUE. 00837 EXIT 00838 END IF 00839 END DO 00840 IF(.NOT.valid_kind) unsetme = .TRUE. 00841 IF(unsetme) bad1(ibend) = 1 00842 END DO 00843 IF(SUM(bad1)/=0) THEN 00844 counter = SIZE(bend_kind_set)-SUM(bad1) 00845 CALL allocate_bend_kind_set(new_bend_kind_set,counter,error) 00846 counter=0 00847 DO ibend=1,SIZE(bend_kind_set) 00848 IF(bad1(ibend)==0) THEN 00849 counter=counter+1 00850 new_bend_kind_set(counter)= bend_kind_set(ibend) 00851 END IF 00852 END DO 00853 counter=0 00854 ALLOCATE(bad2(SIZE(bend_list)),STAT=stat) 00855 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00856 bad2(:)=0 00857 DO ibend=1,SIZE(bend_list) 00858 unsetme = .FALSE. 00859 IF(bend_list(ibend)%bend_kind%id_type==do_ff_undef) unsetme = .TRUE. 00860 IF(bend_list(ibend)%id_type==do_ff_undef) unsetme = .TRUE. 00861 IF(unsetme) bad2(ibend) = 1 00862 END DO 00863 IF(SUM(bad2)/=0) THEN 00864 counter = SIZE(bend_list)-SUM(bad2) 00865 ALLOCATE(new_bend_list(counter),STAT=stat) 00866 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00867 counter=0 00868 DO ibend=1,SIZE(bend_list) 00869 IF(bad2(ibend)==0) THEN 00870 counter=counter+1 00871 new_bend_list(counter) = bend_list(ibend) 00872 newkind = bend_list(ibend)%bend_kind%kind_number 00873 newkind = newkind - SUM(bad1(1:newkind)) 00874 new_bend_list(counter)%bend_kind => new_bend_kind_set(newkind) 00875 END IF 00876 END DO 00877 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00878 nbend=SIZE(new_bend_list),& 00879 bend_kind_set=new_bend_kind_set,& 00880 bend_list=new_bend_list) 00881 DO ibend=1,SIZE(new_bend_kind_set) 00882 new_bend_kind_set(ibend)%kind_number=ibend 00883 END DO 00884 END IF 00885 DEALLOCATE(bad2,STAT=stat) 00886 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00887 CALL deallocate_bend_kind_set(bend_kind_set,error) 00888 DEALLOCATE(bend_list,STAT=stat) 00889 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00890 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New BEND Count: ",& 00891 SIZE(new_bend_list),SIZE(new_bend_kind_set) 00892 IF (iw>0) WRITE(iw,'(3I6)')(new_bend_list(ibend)%a,new_bend_list(ibend)%b,& 00893 new_bend_list(ibend)%c,ibend=1,SIZE(new_bend_list)) 00894 END IF 00895 DEALLOCATE(bad1,STAT=stat) 00896 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00897 END IF 00898 END DO 00899 00900 !----------------------------------------------------------------------------- 00901 !----------------------------------------------------------------------------- 00902 ! 9. Count the number of UNSET Urey-Bradley kinds there are 00903 !----------------------------------------------------------------------------- 00904 DO ikind=1,SIZE(molecule_kind_set) 00905 molecule_kind => molecule_kind_set(ikind) 00906 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00907 nub=nub,& 00908 ub_kind_set=ub_kind_set,& 00909 ub_list=ub_list) 00910 IF(nub>0) THEN 00911 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old UB Count: ",& 00912 SIZE(ub_list),SIZE(ub_kind_set) 00913 IF (iw>0) WRITE(iw,'(3I6)')(ub_list(iub)%a,ub_list(iub)%b,& 00914 ub_list(iub)%c,iub=1,SIZE(ub_list)) 00915 NULLIFY(bad1,bad2) 00916 ALLOCATE(bad1(SIZE(ub_kind_set)),STAT=stat) 00917 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00918 bad1(:)=0 00919 DO iub=1,SIZE(ub_kind_set) 00920 unsetme=.FALSE. 00921 IF(ub_kind_set(iub)%id_type==do_ff_undef) unsetme = .TRUE. 00922 valid_kind=.FALSE. 00923 DO i=1,SIZE(ub_list) 00924 IF(ub_list(i)%id_type/=do_ff_undef.AND.& 00925 ub_list(i)%ub_kind%kind_number==iub) THEN 00926 valid_kind=.TRUE. 00927 EXIT 00928 END IF 00929 END DO 00930 IF(.NOT.valid_kind) unsetme = .TRUE. 00931 IF(unsetme) bad1(iub) = 1 00932 END DO 00933 IF(SUM(bad1)/=0) THEN 00934 counter = SIZE(ub_kind_set)-SUM(bad1) 00935 CALL allocate_ub_kind_set(new_ub_kind_set,counter,error) 00936 counter=0 00937 DO iub=1,SIZE(ub_kind_set) 00938 IF(bad1(iub)==0) THEN 00939 counter=counter+1 00940 new_ub_kind_set(counter)= ub_kind_set(iub) 00941 END IF 00942 END DO 00943 counter=0 00944 ALLOCATE(bad2(SIZE(ub_list)),STAT=stat) 00945 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00946 bad2(:)=0 00947 DO iub=1,SIZE(ub_list) 00948 unsetme = .FALSE. 00949 IF(ub_list(iub)%ub_kind%id_type==do_ff_undef) unsetme = .TRUE. 00950 IF(ub_list(iub)%id_type==do_ff_undef) unsetme = .TRUE. 00951 IF(unsetme) bad2(iub) = 1 00952 END DO 00953 IF(SUM(bad2)/=0) THEN 00954 counter = SIZE(ub_list)-SUM(bad2) 00955 ALLOCATE(new_ub_list(counter),STAT=stat) 00956 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00957 counter=0 00958 DO iub=1,SIZE(ub_list) 00959 IF(bad2(iub)==0) THEN 00960 counter=counter+1 00961 new_ub_list(counter) = ub_list(iub) 00962 newkind = ub_list(iub)%ub_kind%kind_number 00963 newkind = newkind - SUM(bad1(1:newkind)) 00964 new_ub_list(counter)%ub_kind => new_ub_kind_set(newkind) 00965 END IF 00966 END DO 00967 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00968 nub=SIZE(new_ub_list),& 00969 ub_kind_set=new_ub_kind_set,& 00970 ub_list=new_ub_list) 00971 DO iub=1,SIZE(new_ub_kind_set) 00972 new_ub_kind_set(iub)%kind_number=iub 00973 END DO 00974 END IF 00975 DEALLOCATE(bad2,STAT=stat) 00976 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00977 CALL ub_kind_dealloc_ref(ub_kind_set,error=error) 00978 DEALLOCATE(ub_list,STAT=stat) 00979 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00980 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New UB Count: ",& 00981 SIZE(new_ub_list),SIZE(new_ub_kind_set) 00982 IF (iw>0) WRITE(iw,'(3I6)')(new_ub_list(iub)%a,new_ub_list(iub)%b,& 00983 new_ub_list(iub)%c,iub=1,SIZE(new_ub_list)) 00984 END IF 00985 DEALLOCATE(bad1,STAT=stat) 00986 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00987 END IF 00988 END DO 00989 00990 !----------------------------------------------------------------------------- 00991 !----------------------------------------------------------------------------- 00992 ! 10. Count the number of UNSET torsion kinds there are 00993 !----------------------------------------------------------------------------- 00994 DO ikind=1,SIZE(molecule_kind_set) 00995 molecule_kind => molecule_kind_set(ikind) 00996 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00997 ntorsion=ntorsion,& 00998 torsion_kind_set=torsion_kind_set,& 00999 torsion_list=torsion_list) 01000 IF(ntorsion>0) THEN 01001 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old TORSION Count: ",& 01002 SIZE(torsion_list),SIZE(torsion_kind_set) 01003 IF (iw>0) WRITE(iw,'(4I6)')(torsion_list(itorsion)%a,torsion_list(itorsion)%b,& 01004 torsion_list(itorsion)%c,torsion_list(itorsion)%d,itorsion=1,SIZE(torsion_list)) 01005 NULLIFY(bad1,bad2) 01006 ALLOCATE(bad1(SIZE(torsion_kind_set)),STAT=stat) 01007 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01008 bad1(:)=0 01009 DO itorsion=1,SIZE(torsion_kind_set) 01010 unsetme=.FALSE. 01011 IF(torsion_kind_set(itorsion)%id_type==do_ff_undef) unsetme = .TRUE. 01012 valid_kind=.FALSE. 01013 DO i=1,SIZE(torsion_list) 01014 IF(torsion_list(i)%id_type/=do_ff_undef.AND.& 01015 torsion_list(i)%torsion_kind%kind_number==itorsion) THEN 01016 valid_kind=.TRUE. 01017 EXIT 01018 END IF 01019 END DO 01020 IF(.NOT.valid_kind) unsetme = .TRUE. 01021 IF(unsetme) bad1(itorsion) = 1 01022 END DO 01023 IF(SUM(bad1)/=0) THEN 01024 counter = SIZE(torsion_kind_set)-SUM(bad1) 01025 CALL allocate_torsion_kind_set(new_torsion_kind_set,counter,error) 01026 counter=0 01027 DO itorsion=1,SIZE(torsion_kind_set) 01028 IF(bad1(itorsion)==0) THEN 01029 counter=counter+1 01030 new_torsion_kind_set(counter)= torsion_kind_set(itorsion) 01031 i = SIZE(torsion_kind_set(itorsion)%m) 01032 j = SIZE(torsion_kind_set(itorsion)%k) 01033 k = SIZE(torsion_kind_set(itorsion)%phi0) 01034 ALLOCATE(new_torsion_kind_set(counter)%m(i),STAT=stat) 01035 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01036 ALLOCATE(new_torsion_kind_set(counter)%k(i),STAT=stat) 01037 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01038 ALLOCATE(new_torsion_kind_set(counter)%phi0(i),STAT=stat) 01039 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01040 new_torsion_kind_set(counter)%m= torsion_kind_set(itorsion)%m 01041 new_torsion_kind_set(counter)%k= torsion_kind_set(itorsion)%k 01042 new_torsion_kind_set(counter)%phi0= torsion_kind_set(itorsion)%phi0 01043 END IF 01044 END DO 01045 counter=0 01046 ALLOCATE(bad2(SIZE(torsion_list)),STAT=stat) 01047 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01048 bad2(:)=0 01049 DO itorsion=1,SIZE(torsion_list) 01050 unsetme = .FALSE. 01051 IF(torsion_list(itorsion)%torsion_kind%id_type==do_ff_undef) unsetme = .TRUE. 01052 IF(torsion_list(itorsion)%id_type==do_ff_undef) unsetme = .TRUE. 01053 IF(unsetme) bad2(itorsion) = 1 01054 END DO 01055 IF(SUM(bad2)/=0) THEN 01056 counter = SIZE(torsion_list)-SUM(bad2) 01057 ALLOCATE(new_torsion_list(counter),STAT=stat) 01058 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01059 counter=0 01060 DO itorsion=1,SIZE(torsion_list) 01061 IF(bad2(itorsion)==0) THEN 01062 counter=counter+1 01063 new_torsion_list(counter) = torsion_list(itorsion) 01064 newkind = torsion_list(itorsion)%torsion_kind%kind_number 01065 newkind = newkind - SUM(bad1(1:newkind)) 01066 new_torsion_list(counter)%torsion_kind => new_torsion_kind_set(newkind) 01067 END IF 01068 END DO 01069 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01070 ntorsion=SIZE(new_torsion_list),& 01071 torsion_kind_set=new_torsion_kind_set,& 01072 torsion_list=new_torsion_list) 01073 DO itorsion=1,SIZE(new_torsion_kind_set) 01074 new_torsion_kind_set(itorsion)%kind_number=itorsion 01075 END DO 01076 END IF 01077 DEALLOCATE(bad2,STAT=stat) 01078 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01079 DO itorsion=1,SIZE(torsion_kind_set) 01080 CALL torsion_kind_dealloc_ref(torsion_kind_set(itorsion),error=error) 01081 END DO 01082 DEALLOCATE(torsion_kind_set,STAT=stat) 01083 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01084 DEALLOCATE(torsion_list,STAT=stat) 01085 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01086 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New TORSION Count: ",& 01087 SIZE(new_torsion_list),SIZE(new_torsion_kind_set) 01088 IF (iw>0) WRITE(iw,'(4I6)')(new_torsion_list(itorsion)%a,new_torsion_list(itorsion)%b,& 01089 new_torsion_list(itorsion)%c,new_torsion_list(itorsion)%d,itorsion=1,& 01090 SIZE(new_torsion_list)) 01091 END IF 01092 DEALLOCATE(bad1,STAT=stat) 01093 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01094 END IF 01095 END DO 01096 01097 !----------------------------------------------------------------------------- 01098 !----------------------------------------------------------------------------- 01099 ! 11. Count the number of UNSET improper kinds there are 01100 !----------------------------------------------------------------------------- 01101 DO ikind=1,SIZE(molecule_kind_set) 01102 molecule_kind => molecule_kind_set(ikind) 01103 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01104 nimpr=nimpr,& 01105 impr_kind_set=impr_kind_set,& 01106 impr_list=impr_list) 01107 IF(nimpr>0) THEN 01108 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old IMPROPER Count: ",& 01109 SIZE(impr_list),SIZE(impr_kind_set) 01110 NULLIFY(bad1,bad2) 01111 ALLOCATE(bad1(SIZE(impr_kind_set)),STAT=stat) 01112 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01113 bad1(:)=0 01114 DO iimpr=1,SIZE(impr_kind_set) 01115 unsetme=.FALSE. 01116 IF(impr_kind_set(iimpr)%id_type==do_ff_undef) unsetme = .TRUE. 01117 valid_kind=.FALSE. 01118 DO i=1,SIZE(impr_list) 01119 IF(impr_list(i)%id_type/=do_ff_undef.AND.& 01120 impr_list(i)%impr_kind%kind_number==iimpr) THEN 01121 valid_kind=.TRUE. 01122 EXIT 01123 END IF 01124 END DO 01125 IF(.NOT.valid_kind) unsetme = .TRUE. 01126 IF(unsetme) bad1(iimpr) = 1 01127 END DO 01128 IF(SUM(bad1)/=0) THEN 01129 counter = SIZE(impr_kind_set)-SUM(bad1) 01130 CALL allocate_impr_kind_set(new_impr_kind_set,counter,error) 01131 counter=0 01132 DO iimpr=1,SIZE(impr_kind_set) 01133 IF(bad1(iimpr)==0) THEN 01134 counter=counter+1 01135 new_impr_kind_set(counter)= impr_kind_set(iimpr) 01136 END IF 01137 END DO 01138 counter=0 01139 ALLOCATE(bad2(SIZE(impr_list)),STAT=stat) 01140 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01141 bad2(:)=0 01142 DO iimpr=1,SIZE(impr_list) 01143 unsetme = .FALSE. 01144 IF(impr_list(iimpr)%impr_kind%id_type==do_ff_undef) unsetme = .TRUE. 01145 IF(impr_list(iimpr)%id_type==do_ff_undef) unsetme = .TRUE. 01146 IF(unsetme) bad2(iimpr) = 1 01147 END DO 01148 IF(SUM(bad2)/=0) THEN 01149 counter = SIZE(impr_list)-SUM(bad2) 01150 ALLOCATE(new_impr_list(counter),STAT=stat) 01151 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01152 counter=0 01153 DO iimpr=1,SIZE(impr_list) 01154 IF(bad2(iimpr)==0) THEN 01155 counter=counter+1 01156 new_impr_list(counter) = impr_list(iimpr) 01157 newkind = impr_list(iimpr)%impr_kind%kind_number 01158 newkind = newkind - SUM(bad1(1:newkind)) 01159 new_impr_list(counter)%impr_kind => new_impr_kind_set(newkind) 01160 END IF 01161 END DO 01162 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01163 nimpr=SIZE(new_impr_list),& 01164 impr_kind_set=new_impr_kind_set,& 01165 impr_list=new_impr_list) 01166 DO iimpr=1,SIZE(new_impr_kind_set) 01167 new_impr_kind_set(iimpr)%kind_number=iimpr 01168 END DO 01169 END IF 01170 DEALLOCATE(bad2,STAT=stat) 01171 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01172 DO iimpr=1,SIZE(impr_kind_set) 01173 CALL impr_kind_dealloc_ref(impr_kind_set(iimpr),error=error) 01174 END DO 01175 DEALLOCATE(impr_kind_set,STAT=stat) 01176 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01177 DEALLOCATE(impr_list,STAT=stat) 01178 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01179 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New IMPROPER Count: ",& 01180 SIZE(new_impr_list),SIZE(new_impr_kind_set) 01181 END IF 01182 DEALLOCATE(bad1,STAT=stat) 01183 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01184 END IF 01185 END DO 01186 01187 !----------------------------------------------------------------------------- 01188 !----------------------------------------------------------------------------- 01189 ! 11. Count the number of UNSET opbends kinds there are 01190 !----------------------------------------------------------------------------- 01191 DO ikind=1,SIZE(molecule_kind_set) 01192 molecule_kind => molecule_kind_set(ikind) 01193 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01194 nopbend=nopbend,& 01195 opbend_kind_set=opbend_kind_set,& 01196 opbend_list=opbend_list) 01197 IF(nopbend>0) THEN 01198 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") Old OPBEND Count: ",& 01199 SIZE(opbend_list),SIZE(opbend_kind_set) 01200 NULLIFY(bad1,bad2) 01201 ALLOCATE(bad1(SIZE(opbend_kind_set)),STAT=stat) 01202 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01203 bad1(:)=0 01204 DO iopbend=1,SIZE(opbend_kind_set) 01205 unsetme=.FALSE. 01206 IF(opbend_kind_set(iopbend)%id_type==do_ff_undef) unsetme = .TRUE. 01207 valid_kind=.FALSE. 01208 DO i=1,SIZE(opbend_list) 01209 IF(opbend_list(i)%id_type/=do_ff_undef.AND.& 01210 opbend_list(i)%opbend_kind%kind_number==iopbend) THEN 01211 valid_kind=.TRUE. 01212 EXIT 01213 END IF 01214 END DO 01215 IF(.NOT.valid_kind) unsetme = .TRUE. 01216 IF(unsetme) bad1(iopbend) = 1 01217 END DO 01218 IF(SUM(bad1)/=0) THEN 01219 counter = SIZE(opbend_kind_set)-SUM(bad1) 01220 CALL allocate_opbend_kind_set(new_opbend_kind_set,counter,error) 01221 counter=0 01222 DO iopbend=1,SIZE(opbend_kind_set) 01223 IF(bad1(iopbend)==0) THEN 01224 counter=counter+1 01225 new_opbend_kind_set(counter)= opbend_kind_set(iopbend) 01226 END IF 01227 END DO 01228 counter=0 01229 ALLOCATE(bad2(SIZE(opbend_list)),STAT=stat) 01230 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01231 bad2(:)=0 01232 DO iopbend=1,SIZE(opbend_list) 01233 unsetme = .FALSE. 01234 IF(opbend_list(iopbend)%opbend_kind%id_type==do_ff_undef) unsetme = .TRUE. 01235 IF(opbend_list(iopbend)%id_type==do_ff_undef) unsetme = .TRUE. 01236 IF(unsetme) bad2(iopbend) = 1 01237 END DO 01238 IF(SUM(bad2)/=0) THEN 01239 counter = SIZE(opbend_list)-SUM(bad2) 01240 ALLOCATE(new_opbend_list(counter),STAT=stat) 01241 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01242 counter=0 01243 DO iopbend=1,SIZE(opbend_list) 01244 IF(bad2(iopbend)==0) THEN 01245 counter=counter+1 01246 new_opbend_list(counter) = opbend_list(iopbend) 01247 newkind = opbend_list(iopbend)%opbend_kind%kind_number 01248 newkind = newkind - SUM(bad1(1:newkind)) 01249 new_opbend_list(counter)%opbend_kind => new_opbend_kind_set(newkind) 01250 END IF 01251 END DO 01252 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01253 nopbend=SIZE(new_opbend_list),& 01254 opbend_kind_set=new_opbend_kind_set,& 01255 opbend_list=new_opbend_list) 01256 DO iopbend=1,SIZE(new_opbend_kind_set) 01257 new_opbend_kind_set(iopbend)%kind_number=iopbend 01258 END DO 01259 END IF 01260 DEALLOCATE(bad2,STAT=stat) 01261 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01262 DEALLOCATE(opbend_kind_set,STAT=stat) 01263 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01264 DEALLOCATE(opbend_list,STAT=stat) 01265 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01266 IF(iw>0) WRITE(iw,*) " Mol(",ikind,") New OPBEND Count: ",& 01267 SIZE(new_opbend_list),SIZE(new_opbend_kind_set) 01268 END IF 01269 DEALLOCATE(bad1,STAT=stat) 01270 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01271 END IF 01272 END DO 01273 01274 !--------------------------------------------------------------------------- 01275 !--------------------------------------------------------------------------- 01276 ! 12. Count the number of UNSET NONBOND14 kinds there are 01277 !- NEED TO REMOVE EXTRAS HERE - IKUO 01278 !--------------------------------------------------------------------------- 01279 CALL cp_print_key_finished_output(iw,logger,mm_section,& 01280 "PRINT%FF_INFO",error=error) 01281 CALL timestop(handle) 01282 END SUBROUTINE clean_intra_force_kind 01283 01284 ! ***************************************************************************** 01287 SUBROUTINE get_generic_info(gen_section, func_name, xfunction, parameters, values,& 01288 var_values, size_variables, i_rep_sec, input_variables, error) 01289 TYPE(section_vals_type), POINTER :: gen_section 01290 CHARACTER(LEN=*), INTENT(IN) :: func_name 01291 CHARACTER(LEN=default_path_length), 01292 INTENT(OUT) :: xfunction 01293 CHARACTER(LEN=default_string_length), 01294 DIMENSION(:), POINTER :: parameters 01295 REAL(KIND=dp), DIMENSION(:), POINTER :: values 01296 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 01297 POINTER :: var_values 01298 INTEGER, INTENT(IN), OPTIONAL :: size_variables, i_rep_sec 01299 CHARACTER(LEN=*), DIMENSION(:), OPTIONAL :: input_variables 01300 TYPE(cp_error_type), INTENT(inout) :: error 01301 01302 CHARACTER(len=*), PARAMETER :: routineN = 'get_generic_info', 01303 routineP = moduleN//':'//routineN 01304 01305 CHARACTER(LEN=default_string_length), 01306 DIMENSION(:), POINTER :: my_par, my_par_tmp, my_units, 01307 my_units_tmp, my_var 01308 INTEGER :: i, ind, irep, isize, j, 01309 mydim, n_par, n_units, n_val, 01310 nblank, stat 01311 LOGICAL :: check, failure 01312 REAL(KIND=dp), DIMENSION(:), POINTER :: my_val, my_val_tmp 01313 01314 failure = .FALSE. 01315 NULLIFY(my_var, my_par, my_val, my_par_tmp, my_val_tmp) 01316 NULLIFY (my_units) 01317 NULLIFY (my_units_tmp) 01318 IF (ASSOCIATED(parameters)) THEN 01319 DEALLOCATE(parameters,STAT=stat) 01320 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01321 END IF 01322 IF (ASSOCIATED(values)) THEN 01323 DEALLOCATE(values,STAT=stat) 01324 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01325 END IF 01326 irep = 1 01327 IF (PRESENT(i_rep_sec)) irep = i_rep_sec 01328 mydim = 0 01329 CALL section_vals_val_get(gen_section,TRIM(func_name),i_rep_section=irep,c_val=xfunction,error=error) 01330 CALL compress(xfunction, full=.TRUE.) 01331 IF (PRESENT(input_variables)) THEN 01332 ALLOCATE(my_var(SIZE(input_variables)),stat=stat) 01333 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01334 my_var = input_variables 01335 ELSE 01336 CALL section_vals_val_get(gen_section,"VARIABLES",i_rep_section=irep,c_vals=my_var,error=error) 01337 END IF 01338 IF (ASSOCIATED(my_var)) THEN 01339 mydim = SIZE(my_var) 01340 END IF 01341 IF (PRESENT(size_variables)) THEN 01342 CPPrecondition(mydim==size_variables,cp_failure_level,routineP,error,failure) 01343 END IF 01344 ! Check for presence of Parameters 01345 CALL section_vals_val_get(gen_section,"PARAMETERS",i_rep_section=irep,n_rep_val=n_par,error=error) 01346 CALL section_vals_val_get(gen_section,"VALUES",i_rep_section=irep,n_rep_val=n_val,error=error) 01347 check = (n_par>0).EQV.(n_val>0) 01348 CPPrecondition(check,cp_failure_level,routineP,error,failure) 01349 CALL section_vals_val_get(gen_section,"UNITS",i_rep_section=irep,n_rep_val=n_units,error=error) 01350 IF (n_par>0) THEN 01351 ! Parameters 01352 ALLOCATE(my_par(0),stat=stat) 01353 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01354 ALLOCATE(my_val(0),stat=stat) 01355 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01356 DO i = 1, n_par 01357 isize = SIZE(my_par) 01358 CALL section_vals_val_get(gen_section,"PARAMETERS",i_rep_section=irep,i_rep_val=i,c_vals=my_par_tmp,error=error) 01359 nblank = COUNT(my_par_tmp=="") 01360 CALL reallocate(my_par, 1, isize+SIZE(my_par_tmp)-nblank) 01361 ind = 0 01362 DO j = 1, SIZE(my_par_tmp) 01363 IF (my_par_tmp(j)=="") CYCLE 01364 ind = ind + 1 01365 my_par(isize+ind)=my_par_tmp(j) 01366 END DO 01367 END DO 01368 DO i = 1, n_val 01369 isize = SIZE(my_val) 01370 CALL section_vals_val_get(gen_section,"VALUES",i_rep_section=irep,i_rep_val=i,r_vals=my_val_tmp,error=error) 01371 CALL reallocate(my_val,1, isize+SIZE(my_val_tmp)) 01372 my_val(isize+1:isize+SIZE(my_val_tmp))=my_val_tmp 01373 END DO 01374 CPPrecondition(SIZE(my_par)==SIZE(my_val),cp_failure_level,routineP,error,failure) 01375 ! Optionally read the units for each parameter value 01376 ALLOCATE (my_units(0),stat=stat) 01377 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01378 IF (n_units > 0) THEN 01379 DO i=1,n_units 01380 isize = SIZE(my_units) 01381 CALL section_vals_val_get(gen_section,"UNITS",i_rep_section=irep,i_rep_val=i,c_vals=my_units_tmp,error=error) 01382 nblank = COUNT(my_units_tmp == "") 01383 CALL reallocate(my_units,1,isize+SIZE(my_units_tmp)-nblank) 01384 ind = 0 01385 DO j=1,SIZE(my_units_tmp) 01386 IF (my_units_tmp(j) == "") CYCLE 01387 ind = ind + 1 01388 my_units(isize+ind) = my_units_tmp(j) 01389 END DO 01390 END DO 01391 CPPrecondition(SIZE(my_units)==SIZE(my_val),cp_failure_level,routineP,error,failure) 01392 END IF 01393 mydim=mydim+SIZE(my_val) 01394 IF (SIZE(my_val)==0) THEN 01395 DEALLOCATE(my_par,stat=stat) 01396 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01397 DEALLOCATE(my_val,stat=stat) 01398 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01399 DEALLOCATE(my_units,stat=stat) 01400 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01401 END IF 01402 END IF 01403 ! Handle trivial case of a null function defined 01404 ALLOCATE(parameters(mydim),stat=stat) 01405 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01406 ALLOCATE(values(mydim),stat=stat) 01407 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01408 IF (mydim>0) THEN 01409 parameters(1:SIZE(my_var)) = my_var 01410 values(1:SIZE(my_var)) = 0.0_dp 01411 IF (PRESENT(var_values)) THEN 01412 CPPrecondition(SIZE(var_values)==SIZE(my_var),cp_failure_level,routineP,error,failure) 01413 values(1:SIZE(my_var)) = var_values 01414 END IF 01415 IF (ASSOCIATED(my_val)) THEN 01416 DO i=1,SIZE(my_val) 01417 parameters(SIZE(my_var)+i) = my_par(i) 01418 IF (n_units > 0) THEN 01419 values(SIZE(my_var)+i) = cp_unit_to_cp2k(my_val(i),TRIM(ADJUSTL(my_units(i))),error=error) 01420 ELSE 01421 values(SIZE(my_var)+i) = my_val(i) 01422 END IF 01423 END DO 01424 END IF 01425 END IF 01426 IF (ASSOCIATED(my_par)) THEN 01427 DEALLOCATE(my_par,stat=stat) 01428 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01429 END IF 01430 IF (ASSOCIATED(my_val)) THEN 01431 DEALLOCATE(my_val,stat=stat) 01432 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01433 END IF 01434 IF (ASSOCIATED(my_units)) THEN 01435 DEALLOCATE(my_units,stat=stat) 01436 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01437 END IF 01438 IF (PRESENT(input_variables)) THEN 01439 DEALLOCATE(my_var,stat=stat) 01440 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01441 END IF 01442 END SUBROUTINE get_generic_info 01443 01444 END MODULE force_fields_util
1.7.3