|
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 ! ***************************************************************************** 00012 MODULE force_fields_all 00013 USE atomic_kind_types, ONLY: atomic_kind_type,& 00014 get_atomic_kind,& 00015 get_atomic_kind_set,& 00016 set_atomic_kind 00017 USE atoms_input, ONLY: read_shell_coord_input 00018 USE cell_types, ONLY: cell_type 00019 USE cp_linked_list_val, ONLY: cp_sll_val_next,& 00020 cp_sll_val_type 00021 USE damping_dipole_types, ONLY: damping_p_create,& 00022 damping_p_type,& 00023 tang_toennies 00024 USE ewald_environment_types, ONLY: ewald_env_get,& 00025 ewald_env_set,& 00026 ewald_environment_type 00027 USE external_potential_types, ONLY: fist_potential_type,& 00028 get_potential,& 00029 set_potential 00030 USE f77_blas 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, force_field_type, & 00036 gromos_info_type, impr_kind_type, input_info_type, opbend_kind_type, & 00037 torsion_kind_type, ub_kind_type 00038 USE input_constants, ONLY: do_ff_amber,& 00039 do_ff_charmm,& 00040 do_ff_g87,& 00041 do_ff_g96,& 00042 do_ff_undef,& 00043 do_qmmm_none 00044 USE input_cp2k_binary_restarts, ONLY: read_binary_cs_coordinates 00045 USE input_section_types, ONLY: section_vals_get,& 00046 section_vals_get_subs_vals,& 00047 section_vals_list_get,& 00048 section_vals_type,& 00049 section_vals_val_get 00050 USE input_val_types, ONLY: val_get,& 00051 val_type 00052 USE kinds, ONLY: default_path_length,& 00053 default_string_length,& 00054 dp 00055 USE mathconstants, ONLY: sqrthalf 00056 USE memory_utilities, ONLY: reallocate 00057 USE molecule_kind_types, ONLY: & 00058 bend_type, bond_type, get_molecule_kind, impr_type, & 00059 molecule_kind_type, opbend_type, set_molecule_kind, shell_type, & 00060 torsion_type, ub_type 00061 USE molecule_types_new, ONLY: get_molecule,& 00062 molecule_type 00063 USE pair_potential, ONLY: get_nonbond_storage,& 00064 spline_nonbond_control 00065 USE pair_potential_coulomb, ONLY: potential_coulomb 00066 USE pair_potential_types, ONLY: & 00067 ea_type, lj_charmm_type, lj_type, nn_type, nosh_nosh, nosh_sh, & 00068 pair_potential_lj_create, pair_potential_pp_create, & 00069 pair_potential_pp_type, pair_potential_single_add, & 00070 pair_potential_single_clean, pair_potential_single_copy, & 00071 pair_potential_single_type, sh_sh, siepmann_type, tersoff_type 00072 USE particle_types, ONLY: allocate_particle_set,& 00073 particle_type 00074 USE physcon, ONLY: bohr 00075 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 00076 USE qmmm_types, ONLY: qmmm_env_mm_type 00077 USE shell_potential_types, ONLY: shell_create,& 00078 shell_kind_type,& 00079 shell_release 00080 USE splines_types, ONLY: spline_data_p_release,& 00081 spline_data_p_retain,& 00082 spline_data_p_type,& 00083 spline_env_release,& 00084 spline_environment_type 00085 USE string_utilities, ONLY: compress,& 00086 integer_to_string,& 00087 uppercase 00088 USE termination, ONLY: stop_program 00089 USE timings, ONLY: timeset,& 00090 timestop 00091 #include "cp_common_uses.h" 00092 00093 IMPLICIT NONE 00094 00095 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_all' 00096 00097 PRIVATE 00098 PUBLIC :: force_field_unique_bond,& 00099 force_field_unique_bend,& 00100 force_field_unique_ub,& 00101 force_field_unique_tors,& 00102 force_field_unique_impr,& 00103 force_field_unique_opbend,& 00104 force_field_pack_bond,& 00105 force_field_pack_bend,& 00106 force_field_pack_ub,& 00107 force_field_pack_tors,& 00108 force_field_pack_impr,& 00109 force_field_pack_opbend,& 00110 force_field_pack_charge,& 00111 force_field_pack_charges,& 00112 force_field_pack_radius,& 00113 force_field_pack_pol, & 00114 force_field_pack_shell,& 00115 force_field_pack_nonbond14,& 00116 force_field_pack_nonbond,& 00117 force_field_pack_splines,& 00118 force_field_pack_eicut,& 00119 force_field_pack_damp 00120 00121 CONTAINS 00122 00123 ! ***************************************************************************** 00126 SUBROUTINE force_field_unique_bond (particle_set, & 00127 molecule_kind_set, molecule_set, ff_type, error) 00128 00129 TYPE(particle_type), DIMENSION(:), 00130 POINTER :: particle_set 00131 TYPE(molecule_kind_type), DIMENSION(:), 00132 POINTER :: molecule_kind_set 00133 TYPE(molecule_type), DIMENSION(:), 00134 POINTER :: molecule_set 00135 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00136 TYPE(cp_error_type), INTENT(INOUT) :: error 00137 00138 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bond', 00139 routineP = moduleN//':'//routineN 00140 00141 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00142 name_atm_b, name_atm_b2 00143 INTEGER :: atm_a, atm_b, counter, first, 00144 handle2, i, j, k, last, 00145 natom, nbond, stat 00146 INTEGER, DIMENSION(:), POINTER :: molecule_list 00147 INTEGER, POINTER :: map_bond_kind(:) 00148 LOGICAL :: failure, found 00149 TYPE(atomic_kind_type), POINTER :: atomic_kind 00150 TYPE(bond_kind_type), DIMENSION(:), 00151 POINTER :: bond_kind_set 00152 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 00153 TYPE(molecule_kind_type), POINTER :: molecule_kind 00154 TYPE(molecule_type), POINTER :: molecule 00155 00156 failure = .FALSE. 00157 CALL timeset(routineN,handle2) 00158 DO i=1,SIZE(molecule_kind_set) 00159 molecule_kind => molecule_kind_set(i) 00160 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00161 molecule_list=molecule_list,& 00162 natom=natom,& 00163 nbond=nbond,bond_list=bond_list) 00164 molecule=>molecule_set(molecule_list(1)) 00165 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00166 IF(nbond>0) THEN 00167 ALLOCATE(map_bond_kind(nbond),STAT=stat) 00168 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00169 counter=0 00170 IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN 00171 DO j=1,nbond 00172 map_bond_kind(j)=j 00173 END DO 00174 counter=nbond 00175 ELSE 00176 DO j=1,nbond 00177 atm_a = bond_list(j)%a 00178 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00179 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00180 name=name_atm_a) 00181 atm_b = bond_list(j)%b 00182 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00183 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00184 name=name_atm_b) 00185 found = .FALSE. 00186 DO k=1,j-1 00187 atm_a = bond_list(k)%a 00188 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00189 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00190 name=name_atm_a2) 00191 atm_b = bond_list(k)%b 00192 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00193 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00194 name=name_atm_b2) 00195 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00196 ((name_atm_b)==(name_atm_b2))) .OR. & 00197 (((name_atm_a)==(name_atm_b2)) .AND. & 00198 ((name_atm_b)==(name_atm_a2))) ) THEN 00199 found = .TRUE. 00200 map_bond_kind(j) = map_bond_kind(k) 00201 EXIT 00202 END IF 00203 END DO 00204 IF(.NOT.found) THEN 00205 counter=counter+1 00206 map_bond_kind(j) = counter 00207 END IF 00208 END DO 00209 END IF 00210 NULLIFY(bond_kind_set) 00211 CALL allocate_bond_kind_set(bond_kind_set,counter,error) 00212 DO j=1,nbond 00213 bond_list(j)%bond_kind => bond_kind_set(map_bond_kind(j)) 00214 END DO 00215 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00216 bond_kind_set=bond_kind_set,bond_list=bond_list) 00217 DEALLOCATE(map_bond_kind,STAT=stat) 00218 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00219 END IF 00220 END DO 00221 CALL timestop(handle2) 00222 00223 END SUBROUTINE force_field_unique_bond 00224 00225 ! ***************************************************************************** 00228 SUBROUTINE force_field_unique_bend (particle_set, & 00229 molecule_kind_set, molecule_set, ff_type, error) 00230 00231 TYPE(particle_type), DIMENSION(:), 00232 POINTER :: particle_set 00233 TYPE(molecule_kind_type), DIMENSION(:), 00234 POINTER :: molecule_kind_set 00235 TYPE(molecule_type), DIMENSION(:), 00236 POINTER :: molecule_set 00237 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00238 TYPE(cp_error_type), INTENT(INOUT) :: error 00239 00240 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bend', 00241 routineP = moduleN//':'//routineN 00242 00243 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00244 name_atm_b, name_atm_b2, 00245 name_atm_c, name_atm_c2 00246 INTEGER :: atm_a, atm_b, atm_c, counter, 00247 first, handle2, i, j, k, 00248 last, natom, nbend, stat 00249 INTEGER, DIMENSION(:), POINTER :: molecule_list 00250 INTEGER, POINTER :: map_bend_kind(:) 00251 LOGICAL :: failure, found 00252 TYPE(atomic_kind_type), POINTER :: atomic_kind 00253 TYPE(bend_kind_type), DIMENSION(:), 00254 POINTER :: bend_kind_set 00255 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 00256 TYPE(molecule_kind_type), POINTER :: molecule_kind 00257 TYPE(molecule_type), POINTER :: molecule 00258 00259 failure = .FALSE. 00260 00261 CALL timeset(routineN,handle2) 00262 DO i=1,SIZE(molecule_kind_set) 00263 molecule_kind => molecule_kind_set(i) 00264 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00265 molecule_list=molecule_list,& 00266 natom=natom,& 00267 nbend=nbend,bend_list=bend_list) 00268 molecule=>molecule_set(molecule_list(1)) 00269 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00270 IF(nbend>0) THEN 00271 ALLOCATE(map_bend_kind(nbend),STAT=stat) 00272 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00273 counter=0 00274 IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN 00275 DO j=1,nbend 00276 map_bend_kind(j)=j 00277 END DO 00278 counter=nbend 00279 ELSE 00280 DO j=1,nbend 00281 atm_a = bend_list(j)%a 00282 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00283 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00284 name=name_atm_a) 00285 atm_b = bend_list(j)%b 00286 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00287 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00288 name=name_atm_b) 00289 atm_c = bend_list(j)%c 00290 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00291 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00292 name=name_atm_c) 00293 found = .FALSE. 00294 DO k=1,j-1 00295 atm_a = bend_list(k)%a 00296 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00297 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00298 name=name_atm_a2) 00299 atm_b = bend_list(k)%b 00300 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00301 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00302 name=name_atm_b2) 00303 atm_c = bend_list(k)%c 00304 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00305 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00306 name=name_atm_c2) 00307 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00308 ((name_atm_b)==(name_atm_b2)) .AND. & 00309 ((name_atm_c)==(name_atm_c2))) .OR. & 00310 (((name_atm_a)==(name_atm_c2)) .AND. & 00311 ((name_atm_b)==(name_atm_b2)) .AND. & 00312 ((name_atm_c)==(name_atm_a2))) ) THEN 00313 found = .TRUE. 00314 map_bend_kind(j) = map_bend_kind(k) 00315 EXIT 00316 END IF 00317 END DO 00318 IF(.NOT.found) THEN 00319 counter=counter+1 00320 map_bend_kind(j) = counter 00321 END IF 00322 END DO 00323 END IF 00324 NULLIFY(bend_kind_set) 00325 CALL allocate_bend_kind_set(bend_kind_set,counter,error) 00326 DO j=1,nbend 00327 bend_list(j)%bend_kind => bend_kind_set(map_bend_kind(j)) 00328 END DO 00329 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00330 bend_kind_set=bend_kind_set,bend_list=bend_list) 00331 DEALLOCATE(map_bend_kind,STAT=stat) 00332 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00333 END IF 00334 END DO 00335 CALL timestop(handle2) 00336 00337 END SUBROUTINE force_field_unique_bend 00338 00339 ! ***************************************************************************** 00342 SUBROUTINE force_field_unique_ub(particle_set, & 00343 molecule_kind_set, molecule_set, error) 00344 00345 TYPE(particle_type), DIMENSION(:), 00346 POINTER :: particle_set 00347 TYPE(molecule_kind_type), DIMENSION(:), 00348 POINTER :: molecule_kind_set 00349 TYPE(molecule_type), DIMENSION(:), 00350 POINTER :: molecule_set 00351 TYPE(cp_error_type), INTENT(INOUT) :: error 00352 00353 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_ub', 00354 routineP = moduleN//':'//routineN 00355 00356 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00357 name_atm_b, name_atm_b2, 00358 name_atm_c, name_atm_c2 00359 INTEGER :: atm_a, atm_b, atm_c, counter, 00360 first, handle2, i, j, k, 00361 last, natom, nub, stat 00362 INTEGER, DIMENSION(:), POINTER :: molecule_list 00363 INTEGER, POINTER :: map_ub_kind(:) 00364 LOGICAL :: failure, found 00365 TYPE(atomic_kind_type), POINTER :: atomic_kind 00366 TYPE(molecule_kind_type), POINTER :: molecule_kind 00367 TYPE(molecule_type), POINTER :: molecule 00368 TYPE(ub_kind_type), DIMENSION(:), 00369 POINTER :: ub_kind_set 00370 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 00371 00372 failure = .FALSE. 00373 00374 CALL timeset(routineN,handle2) 00375 DO i=1,SIZE(molecule_kind_set) 00376 molecule_kind => molecule_kind_set(i) 00377 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00378 molecule_list=molecule_list,& 00379 natom=natom,& 00380 nub=nub,ub_list=ub_list) 00381 molecule=>molecule_set(molecule_list(1)) 00382 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00383 IF(nub>0) THEN 00384 ALLOCATE(map_ub_kind(nub),STAT=stat) 00385 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00386 counter=0 00387 DO j=1,nub 00388 atm_a = ub_list(j)%a 00389 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00390 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00391 name=name_atm_a) 00392 atm_b = ub_list(j)%b 00393 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00394 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00395 name=name_atm_b) 00396 atm_c = ub_list(j)%c 00397 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00398 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00399 name=name_atm_c) 00400 found = .FALSE. 00401 DO k=1,j-1 00402 atm_a = ub_list(k)%a 00403 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00404 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00405 name=name_atm_a2) 00406 atm_b = ub_list(k)%b 00407 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00408 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00409 name=name_atm_b2) 00410 atm_c = ub_list(k)%c 00411 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00412 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00413 name=name_atm_c2) 00414 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00415 ((name_atm_b)==(name_atm_b2)) .AND. & 00416 ((name_atm_c)==(name_atm_c2))) .OR. & 00417 (((name_atm_a)==(name_atm_c2)) .AND. & 00418 ((name_atm_b)==(name_atm_b2)) .AND. & 00419 ((name_atm_c)==(name_atm_a2))) ) THEN 00420 found = .TRUE. 00421 map_ub_kind(j) = map_ub_kind(k) 00422 EXIT 00423 END IF 00424 END DO 00425 IF(.NOT.found) THEN 00426 counter=counter+1 00427 map_ub_kind(j) = counter 00428 END IF 00429 END DO 00430 CALL allocate_ub_kind_set(ub_kind_set,counter,error) 00431 DO j=1,nub 00432 ub_list(j)%ub_kind => ub_kind_set(map_ub_kind(j)) 00433 END DO 00434 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00435 ub_kind_set=ub_kind_set,ub_list=ub_list) 00436 DEALLOCATE(map_ub_kind,STAT=stat) 00437 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00438 END IF 00439 END DO 00440 CALL timestop(handle2) 00441 00442 END SUBROUTINE force_field_unique_ub 00443 00444 ! ***************************************************************************** 00447 SUBROUTINE force_field_unique_tors(particle_set, & 00448 molecule_kind_set, molecule_set, ff_type, error) 00449 00450 TYPE(particle_type), DIMENSION(:), 00451 POINTER :: particle_set 00452 TYPE(molecule_kind_type), DIMENSION(:), 00453 POINTER :: molecule_kind_set 00454 TYPE(molecule_type), DIMENSION(:), 00455 POINTER :: molecule_set 00456 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00457 TYPE(cp_error_type), INTENT(INOUT) :: error 00458 00459 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_tors', 00460 routineP = moduleN//':'//routineN 00461 00462 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00463 name_atm_b, name_atm_b2, name_atm_c, name_atm_c2, name_atm_d, 00464 name_atm_d2 00465 INTEGER :: atm_a, atm_b, atm_c, atm_d, 00466 counter, first, handle2, i, 00467 j, k, last, natom, ntorsion, 00468 stat 00469 INTEGER, DIMENSION(:), POINTER :: molecule_list 00470 INTEGER, POINTER :: map_torsion_kind(:) 00471 LOGICAL :: failure, found 00472 TYPE(atomic_kind_type), POINTER :: atomic_kind 00473 TYPE(molecule_kind_type), POINTER :: molecule_kind 00474 TYPE(molecule_type), POINTER :: molecule 00475 TYPE(torsion_kind_type), DIMENSION(:), 00476 POINTER :: torsion_kind_set 00477 TYPE(torsion_type), DIMENSION(:), 00478 POINTER :: torsion_list 00479 00480 failure = .FALSE. 00481 CALL timeset(routineN,handle2) 00482 DO i=1,SIZE(molecule_kind_set) 00483 molecule_kind => molecule_kind_set(i) 00484 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00485 molecule_list=molecule_list,& 00486 natom=natom,& 00487 ntorsion=ntorsion,torsion_list=torsion_list) 00488 molecule=>molecule_set(molecule_list(1)) 00489 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00490 IF(ntorsion>0) THEN 00491 ALLOCATE(map_torsion_kind(ntorsion),STAT=stat) 00492 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00493 counter=0 00494 IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN 00495 DO j=1,ntorsion 00496 map_torsion_kind(j)=j 00497 END DO 00498 counter=ntorsion 00499 ELSE 00500 DO j=1,ntorsion 00501 atm_a = torsion_list(j)%a 00502 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00503 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00504 name=name_atm_a) 00505 atm_b = torsion_list(j)%b 00506 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00507 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00508 name=name_atm_b) 00509 atm_c = torsion_list(j)%c 00510 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00511 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00512 name=name_atm_c) 00513 atm_d = torsion_list(j)%d 00514 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00515 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00516 name=name_atm_d) 00517 found = .FALSE. 00518 DO k=1,j-1 00519 atm_a = torsion_list(k)%a 00520 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00521 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00522 name=name_atm_a2) 00523 atm_b = torsion_list(k)%b 00524 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00525 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00526 name=name_atm_b2) 00527 atm_c = torsion_list(k)%c 00528 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00529 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00530 name=name_atm_c2) 00531 atm_d = torsion_list(k)%d 00532 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00533 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00534 name=name_atm_d2) 00535 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00536 ((name_atm_b)==(name_atm_b2)) .AND. & 00537 ((name_atm_c)==(name_atm_c2)) .AND. & 00538 ((name_atm_d)==(name_atm_d2))) .OR. & 00539 (((name_atm_a)==(name_atm_d2)) .AND. & 00540 ((name_atm_b)==(name_atm_c2)) .AND. & 00541 ((name_atm_c)==(name_atm_b2)) .AND. & 00542 ((name_atm_d)==(name_atm_a2))) ) THEN 00543 found = .TRUE. 00544 map_torsion_kind(j) = map_torsion_kind(k) 00545 EXIT 00546 END IF 00547 END DO 00548 IF(.NOT.found) THEN 00549 counter=counter+1 00550 map_torsion_kind(j) = counter 00551 END IF 00552 END DO 00553 END IF 00554 NULLIFY(torsion_kind_set) 00555 CALL allocate_torsion_kind_set(torsion_kind_set,counter,error) 00556 DO j=1,ntorsion 00557 torsion_list(j)%torsion_kind => torsion_kind_set(map_torsion_kind(j)) 00558 END DO 00559 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00560 torsion_kind_set=torsion_kind_set,torsion_list=torsion_list) 00561 DEALLOCATE(map_torsion_kind,STAT=stat) 00562 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00563 END IF 00564 END DO 00565 CALL timestop(handle2) 00566 00567 END SUBROUTINE force_field_unique_tors 00568 00569 ! ***************************************************************************** 00572 SUBROUTINE force_field_unique_impr (particle_set, & 00573 molecule_kind_set, molecule_set, ff_type, error) 00574 00575 TYPE(particle_type), DIMENSION(:), 00576 POINTER :: particle_set 00577 TYPE(molecule_kind_type), DIMENSION(:), 00578 POINTER :: molecule_kind_set 00579 TYPE(molecule_type), DIMENSION(:), 00580 POINTER :: molecule_set 00581 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00582 TYPE(cp_error_type), INTENT(INOUT) :: error 00583 00584 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_impr', 00585 routineP = moduleN//':'//routineN 00586 00587 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00588 name_atm_b, name_atm_b2, name_atm_c, name_atm_c2, name_atm_d, 00589 name_atm_d2 00590 INTEGER :: atm_a, atm_b, atm_c, atm_d, 00591 counter, first, handle2, i, 00592 j, k, last, natom, nimpr, stat 00593 INTEGER, DIMENSION(:), POINTER :: molecule_list 00594 INTEGER, POINTER :: map_impr_kind(:) 00595 LOGICAL :: failure, found 00596 TYPE(atomic_kind_type), POINTER :: atomic_kind 00597 TYPE(impr_kind_type), DIMENSION(:), 00598 POINTER :: impr_kind_set 00599 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 00600 TYPE(molecule_kind_type), POINTER :: molecule_kind 00601 TYPE(molecule_type), POINTER :: molecule 00602 00603 failure = .FALSE. 00604 00605 CALL timeset(routineN,handle2) 00606 DO i=1,SIZE(molecule_kind_set) 00607 molecule_kind => molecule_kind_set(i) 00608 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00609 molecule_list=molecule_list,& 00610 natom=natom,& 00611 nimpr=nimpr,impr_list=impr_list) 00612 molecule=>molecule_set(molecule_list(1)) 00613 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00614 IF(nimpr>0) THEN 00615 ALLOCATE(map_impr_kind(nimpr),STAT=stat) 00616 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00617 counter=0 00618 IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN 00619 DO j=1,nimpr 00620 map_impr_kind(j)=j 00621 END DO 00622 counter=nimpr 00623 ELSE 00624 DO j=1,nimpr 00625 atm_a = impr_list(j)%a 00626 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00627 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00628 name=name_atm_a) 00629 atm_b = impr_list(j)%b 00630 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00631 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00632 name=name_atm_b) 00633 atm_c = impr_list(j)%c 00634 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00635 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00636 name=name_atm_c) 00637 atm_d = impr_list(j)%d 00638 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00639 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00640 name=name_atm_d) 00641 found = .FALSE. 00642 DO k=1,j-1 00643 atm_a = impr_list(k)%a 00644 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00645 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00646 name=name_atm_a2) 00647 atm_b = impr_list(k)%b 00648 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00649 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00650 name=name_atm_b2) 00651 atm_c = impr_list(k)%c 00652 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00653 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00654 name=name_atm_c2) 00655 atm_d = impr_list(k)%d 00656 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00657 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00658 name=name_atm_d2) 00659 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00660 ((name_atm_b)==(name_atm_b2)) .AND. & 00661 ((name_atm_c)==(name_atm_c2)) .AND. & 00662 ((name_atm_d)==(name_atm_d2))) .OR. & 00663 (((name_atm_a)==(name_atm_a2)) .AND. & 00664 ((name_atm_b)==(name_atm_b2)) .AND. & 00665 ((name_atm_c)==(name_atm_d2)) .AND. & 00666 ((name_atm_d)==(name_atm_c2))) ) THEN 00667 found = .TRUE. 00668 map_impr_kind(j) = map_impr_kind(k) 00669 EXIT 00670 END IF 00671 END DO 00672 IF(.NOT.found) THEN 00673 counter=counter+1 00674 map_impr_kind(j) = counter 00675 END IF 00676 END DO 00677 END IF 00678 NULLIFY(impr_kind_set) 00679 CALL allocate_impr_kind_set(impr_kind_set,counter,error) 00680 DO j=1,nimpr 00681 impr_list(j)%impr_kind => impr_kind_set(map_impr_kind(j)) 00682 END DO 00683 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00684 impr_kind_set=impr_kind_set,impr_list=impr_list) 00685 DEALLOCATE(map_impr_kind,STAT=stat) 00686 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00687 END IF 00688 END DO 00689 CALL timestop(handle2) 00690 00691 END SUBROUTINE force_field_unique_impr 00692 00693 ! ***************************************************************************** 00698 SUBROUTINE force_field_unique_opbend (particle_set, & 00699 molecule_kind_set, molecule_set, ff_type, error) 00700 00701 TYPE(particle_type), DIMENSION(:), 00702 POINTER :: particle_set 00703 TYPE(molecule_kind_type), DIMENSION(:), 00704 POINTER :: molecule_kind_set 00705 TYPE(molecule_type), DIMENSION(:), 00706 POINTER :: molecule_set 00707 TYPE(force_field_type), INTENT(INOUT) :: ff_type 00708 TYPE(cp_error_type), INTENT(INOUT) :: error 00709 00710 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_opbend', 00711 routineP = moduleN//':'//routineN 00712 00713 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, 00714 name_atm_b, name_atm_b2, name_atm_c, name_atm_c2, name_atm_d, 00715 name_atm_d2 00716 INTEGER :: atm_a, atm_b, atm_c, atm_d, 00717 counter, first, handle2, i, 00718 j, k, last, natom, nopbend, 00719 stat 00720 INTEGER, DIMENSION(:), POINTER :: molecule_list 00721 INTEGER, POINTER :: map_opbend_kind(:) 00722 LOGICAL :: failure, found 00723 TYPE(atomic_kind_type), POINTER :: atomic_kind 00724 TYPE(molecule_kind_type), POINTER :: molecule_kind 00725 TYPE(molecule_type), POINTER :: molecule 00726 TYPE(opbend_kind_type), DIMENSION(:), 00727 POINTER :: opbend_kind_set 00728 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 00729 00730 failure = .FALSE. 00731 00732 CALL timeset(routineN,handle2) 00733 DO i=1,SIZE(molecule_kind_set) 00734 molecule_kind => molecule_kind_set(i) 00735 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00736 molecule_list=molecule_list,& 00737 natom=natom,& 00738 nopbend=nopbend,opbend_list=opbend_list) 00739 molecule=>molecule_set(molecule_list(1)) 00740 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00741 IF(nopbend>0) THEN 00742 ALLOCATE(map_opbend_kind(nopbend),STAT=stat) 00743 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00744 counter=0 00745 IF((ff_type%ff_type==do_ff_g96).OR.(ff_type%ff_type==do_ff_g87)) THEN 00746 DO j=1,nopbend 00747 map_opbend_kind(j)=j 00748 END DO 00749 counter=nopbend 00750 ELSE 00751 DO j=1,nopbend 00752 atm_a = opbend_list(j)%a 00753 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00754 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00755 name=name_atm_a) 00756 atm_b = opbend_list(j)%b 00757 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00758 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00759 name=name_atm_b) 00760 atm_c = opbend_list(j)%c 00761 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00762 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00763 name=name_atm_c) 00764 atm_d = opbend_list(j)%d 00765 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00766 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00767 name=name_atm_d) 00768 found = .FALSE. 00769 DO k=1,j-1 00770 atm_a = opbend_list(k)%a 00771 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00772 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00773 name=name_atm_a2) 00774 atm_b = opbend_list(k)%b 00775 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00776 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00777 name=name_atm_b2) 00778 atm_c = opbend_list(k)%c 00779 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 00780 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00781 name=name_atm_c2) 00782 atm_d = opbend_list(k)%d 00783 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 00784 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00785 name=name_atm_d2) 00786 IF( (((name_atm_a)==(name_atm_a2)) .AND. & 00787 ((name_atm_b)==(name_atm_b2)) .AND. & 00788 ((name_atm_c)==(name_atm_c2)) .AND. & 00789 ((name_atm_d)==(name_atm_d2))) .OR. & 00790 (((name_atm_a)==(name_atm_a2)) .AND. & 00791 ((name_atm_b)==(name_atm_c2)) .AND. & 00792 ((name_atm_c)==(name_atm_b2)) .AND. & 00793 ((name_atm_d)==(name_atm_d2))) ) THEN 00794 found = .TRUE. 00795 map_opbend_kind(j) = map_opbend_kind(k) 00796 EXIT 00797 END IF 00798 END DO 00799 IF(.NOT.found) THEN 00800 counter=counter+1 00801 map_opbend_kind(j) = counter 00802 END IF 00803 END DO 00804 END IF 00805 NULLIFY(opbend_kind_set) 00806 CALL allocate_opbend_kind_set(opbend_kind_set,counter,error) 00807 DO j=1,nopbend 00808 opbend_list(j)%opbend_kind => opbend_kind_set(map_opbend_kind(j)) 00809 END DO 00810 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00811 opbend_kind_set=opbend_kind_set,opbend_list=opbend_list) 00812 DEALLOCATE(map_opbend_kind,STAT=stat) 00813 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00814 END IF 00815 END DO 00816 CALL timestop(handle2) 00817 00818 END SUBROUTINE force_field_unique_opbend 00819 00820 ! ***************************************************************************** 00823 SUBROUTINE force_field_pack_bond (particle_set, molecule_kind_set, molecule_set, & 00824 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) 00825 00826 TYPE(particle_type), DIMENSION(:), 00827 POINTER :: particle_set 00828 TYPE(molecule_kind_type), DIMENSION(:), 00829 POINTER :: molecule_kind_set 00830 TYPE(molecule_type), DIMENSION(:), 00831 POINTER :: molecule_set 00832 LOGICAL :: fatal 00833 CHARACTER(LEN=default_string_length), 00834 DIMENSION(:), POINTER :: Ainfo 00835 TYPE(charmm_info_type), POINTER :: chm_info 00836 TYPE(input_info_type), POINTER :: inp_info 00837 TYPE(gromos_info_type), POINTER :: gro_info 00838 TYPE(amber_info_type), POINTER :: amb_info 00839 TYPE(cp_error_type), INTENT(INOUT) :: error 00840 00841 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bond', 00842 routineP = moduleN//':'//routineN 00843 00844 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b 00845 INTEGER :: atm_a, atm_b, first, handle2, 00846 i, itype, j, k, last, natom, 00847 nbond 00848 INTEGER, DIMENSION(:), POINTER :: molecule_list 00849 LOGICAL :: failure, found, only_qm 00850 TYPE(atomic_kind_type), POINTER :: atomic_kind 00851 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 00852 TYPE(molecule_kind_type), POINTER :: molecule_kind 00853 TYPE(molecule_type), POINTER :: molecule 00854 00855 failure = .FALSE. 00856 CALL timeset(routineN,handle2) 00857 00858 DO i=1,SIZE(molecule_kind_set) 00859 molecule_kind => molecule_kind_set(i) 00860 CALL get_molecule_kind(molecule_kind=molecule_kind,& 00861 molecule_list=molecule_list,& 00862 natom=natom,& 00863 nbond=nbond,bond_list=bond_list) 00864 molecule=>molecule_set(molecule_list(1)) 00865 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 00866 DO j=1,nbond 00867 atm_a = bond_list(j)%a 00868 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 00869 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00870 name=name_atm_a) 00871 atm_b = bond_list(j)%b 00872 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 00873 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00874 name=name_atm_b) 00875 found = .FALSE. 00876 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 00877 CALL uppercase(name_atm_a) 00878 CALL uppercase(name_atm_b) 00879 00880 ! loop over params from GROMOS 00881 IF(ASSOCIATED(gro_info%bond_k)) THEN 00882 k=SIZE(gro_info%bond_k) 00883 itype = bond_list(j)%itype 00884 IF(itype<=k) THEN 00885 bond_list(j)%bond_kind%k(1) = gro_info%bond_k(itype) 00886 bond_list(j)%bond_kind%r0 = gro_info%bond_r0(itype) 00887 ELSE 00888 itype=itype-k 00889 bond_list(j)%bond_kind%k(1) = gro_info%solvent_k(itype) 00890 bond_list(j)%bond_kind%r0 = gro_info%solvent_r0(itype) 00891 END IF 00892 bond_list(j)%bond_kind%id_type = gro_info%ff_gromos_type 00893 bond_list(j)%id_type = gro_info%ff_gromos_type 00894 found = .TRUE. 00895 END IF 00896 00897 ! loop over params from CHARMM 00898 IF(ASSOCIATED(chm_info%bond_a)) THEN 00899 DO k=1,SIZE(chm_info%bond_a) 00900 IF( (((chm_info%bond_a(k))==(name_atm_a)) .AND. & 00901 ((chm_info%bond_b(k))==(name_atm_b))) .OR. & 00902 (((chm_info%bond_a(k))==(name_atm_b)) .AND. & 00903 ((chm_info%bond_b(k))==(name_atm_a))) ) THEN 00904 bond_list(j)%bond_kind%id_type = do_ff_charmm 00905 bond_list(j)%bond_kind%k(1) = chm_info%bond_k(k) 00906 bond_list(j)%bond_kind%r0 = chm_info%bond_r0(k) 00907 CALL issue_duplications(found,routineP,"Bond", name_atm_a, name_atm_b,& 00908 error=error) 00909 found = .TRUE. 00910 EXIT 00911 END IF 00912 END DO 00913 END IF 00914 00915 ! loop over params from AMBER 00916 IF(ASSOCIATED(amb_info%bond_a)) THEN 00917 DO k=1,SIZE(amb_info%bond_a) 00918 IF( (((amb_info%bond_a(k))==(name_atm_a)) .AND. & 00919 ((amb_info%bond_b(k))==(name_atm_b))) .OR. & 00920 (((amb_info%bond_a(k))==(name_atm_b)) .AND. & 00921 ((amb_info%bond_b(k))==(name_atm_a))) ) THEN 00922 bond_list(j)%bond_kind%id_type = do_ff_amber 00923 bond_list(j)%bond_kind%k(1) = amb_info%bond_k(k) 00924 bond_list(j)%bond_kind%r0 = amb_info%bond_r0(k) 00925 CALL issue_duplications(found,routineP,"Bond", name_atm_a, name_atm_b,& 00926 error=error) 00927 found = .TRUE. 00928 EXIT 00929 END IF 00930 END DO 00931 END IF 00932 00933 ! always have the input param last to overwrite all the other ones 00934 IF(ASSOCIATED(inp_info%bond_a)) THEN 00935 DO k=1,SIZE(inp_info%bond_a) 00936 IF( (((inp_info%bond_a(k))==(name_atm_a)) .AND. & 00937 ((inp_info%bond_b(k))==(name_atm_b))) .OR. & 00938 (((inp_info%bond_a(k))==(name_atm_b)) .AND. & 00939 ((inp_info%bond_b(k))==(name_atm_a))) ) THEN 00940 bond_list(j)%bond_kind%id_type = inp_info%bond_kind(k) 00941 bond_list(j)%bond_kind%k(:) = inp_info%bond_k(:,k) 00942 bond_list(j)%bond_kind%r0 = inp_info%bond_r0(k) 00943 bond_list(j)%bond_kind%cs = inp_info%bond_cs(k) 00944 CALL issue_duplications(found,routineP,"Bond", name_atm_a, name_atm_b,& 00945 error=error) 00946 found = .TRUE. 00947 EXIT 00948 END IF 00949 END DO 00950 END IF 00951 00952 IF(.NOT.found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 00953 atm2=TRIM(name_atm_b),& 00954 fatal=fatal,& 00955 type_name="Bond",& 00956 array=Ainfo,& 00957 error=error) 00958 ! QM/MM modifications 00959 IF (only_qm) THEN 00960 bond_list(j)%id_type = do_ff_undef 00961 bond_list(j)%bond_kind%id_type = do_ff_undef 00962 END IF 00963 END DO 00964 CALL set_molecule_kind(molecule_kind=molecule_kind,& 00965 bond_list=bond_list) 00966 END DO 00967 CALL timestop(handle2) 00968 00969 END SUBROUTINE force_field_pack_bond 00970 00971 ! ***************************************************************************** 00974 SUBROUTINE force_field_pack_bend (particle_set, molecule_kind_set, molecule_set, & 00975 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info, error) 00976 00977 TYPE(particle_type), DIMENSION(:), 00978 POINTER :: particle_set 00979 TYPE(molecule_kind_type), DIMENSION(:), 00980 POINTER :: molecule_kind_set 00981 TYPE(molecule_type), DIMENSION(:), 00982 POINTER :: molecule_set 00983 LOGICAL :: fatal 00984 CHARACTER(LEN=default_string_length), 00985 DIMENSION(:), POINTER :: Ainfo 00986 TYPE(charmm_info_type), POINTER :: chm_info 00987 TYPE(input_info_type), POINTER :: inp_info 00988 TYPE(gromos_info_type), POINTER :: gro_info 00989 TYPE(amber_info_type), POINTER :: amb_info 00990 TYPE(cp_error_type), INTENT(INOUT) :: error 00991 00992 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bend', 00993 routineP = moduleN//':'//routineN 00994 00995 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, 00996 name_atm_c 00997 INTEGER :: atm_a, atm_b, atm_c, first, 00998 handle2, i, itype, j, k, 00999 last, natom, nbend 01000 INTEGER, DIMENSION(:), POINTER :: molecule_list 01001 LOGICAL :: failure, found, only_qm 01002 TYPE(atomic_kind_type), POINTER :: atomic_kind 01003 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 01004 TYPE(molecule_kind_type), POINTER :: molecule_kind 01005 TYPE(molecule_type), POINTER :: molecule 01006 01007 CALL timeset(routineN,handle2) 01008 failure = .FALSE. 01009 01010 DO i=1,SIZE(molecule_kind_set) 01011 molecule_kind => molecule_kind_set(i) 01012 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01013 molecule_list=molecule_list,& 01014 natom=natom,& 01015 nbend=nbend,bend_list=bend_list) 01016 molecule=>molecule_set(molecule_list(1)) 01017 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 01018 DO j=1,nbend 01019 atm_a = bend_list(j)%a 01020 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 01021 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01022 name=name_atm_a) 01023 atm_b = bend_list(j)%b 01024 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 01025 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01026 name=name_atm_b) 01027 atm_c = bend_list(j)%c 01028 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 01029 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01030 name=name_atm_c) 01031 found = .FALSE. 01032 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c) 01033 CALL uppercase(name_atm_a) 01034 CALL uppercase(name_atm_b) 01035 CALL uppercase(name_atm_c) 01036 01037 ! loop over params from GROMOS 01038 IF(ASSOCIATED(gro_info%bend_k)) THEN 01039 k=SIZE(gro_info%bend_k) 01040 itype = bend_list(j)%itype 01041 IF(itype>0) THEN 01042 bend_list(j)%bend_kind%k = gro_info%bend_k(itype) 01043 bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype) 01044 ELSE 01045 bend_list(j)%bend_kind%k = gro_info%bend_k(itype/k) 01046 bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype/k) 01047 END IF 01048 bend_list(j)%bend_kind%id_type = gro_info%ff_gromos_type 01049 bend_list(j)%id_type = gro_info%ff_gromos_type 01050 found = .TRUE. 01051 END IF 01052 01053 ! loop over params from CHARMM 01054 IF(ASSOCIATED(chm_info%bend_a)) THEN 01055 DO k=1,SIZE(chm_info%bend_a) 01056 IF( (((chm_info%bend_a(k))==(name_atm_a)) .AND. & 01057 ((chm_info%bend_b(k))==(name_atm_b)) .AND. & 01058 ((chm_info%bend_c(k))==(name_atm_c))) .OR. & 01059 (((chm_info%bend_a(k))==(name_atm_c)) .AND. & 01060 ((chm_info%bend_b(k))==(name_atm_b)) .AND. & 01061 ((chm_info%bend_c(k))==(name_atm_a))) ) THEN 01062 bend_list(j)%bend_kind%id_type = do_ff_charmm 01063 bend_list(j)%bend_kind%k = chm_info%bend_k(k) 01064 bend_list(j)%bend_kind%theta0 = chm_info%bend_theta0(k) 01065 CALL issue_duplications(found,routineP, "Bend", name_atm_a, name_atm_b,& 01066 name_atm_c, error=error) 01067 found = .TRUE. 01068 EXIT 01069 END IF 01070 END DO 01071 END IF 01072 01073 ! loop over params from AMBER 01074 IF(ASSOCIATED(amb_info%bend_a)) THEN 01075 DO k=1,SIZE(amb_info%bend_a) 01076 IF( (((amb_info%bend_a(k))==(name_atm_a)) .AND. & 01077 ((amb_info%bend_b(k))==(name_atm_b)) .AND. & 01078 ((amb_info%bend_c(k))==(name_atm_c))) .OR. & 01079 (((amb_info%bend_a(k))==(name_atm_c)) .AND. & 01080 ((amb_info%bend_b(k))==(name_atm_b)) .AND. & 01081 ((amb_info%bend_c(k))==(name_atm_a))) ) THEN 01082 bend_list(j)%bend_kind%id_type = do_ff_amber 01083 bend_list(j)%bend_kind%k = amb_info%bend_k(k) 01084 bend_list(j)%bend_kind%theta0 = amb_info%bend_theta0(k) 01085 CALL issue_duplications(found,routineP, "Bend", name_atm_a, name_atm_b,& 01086 name_atm_c, error=error) 01087 found = .TRUE. 01088 EXIT 01089 END IF 01090 END DO 01091 END IF 01092 01093 ! always have the input param last to overwrite all the other ones 01094 IF(ASSOCIATED(inp_info%bend_a)) THEN 01095 DO k=1,SIZE(inp_info%bend_a) 01096 IF( (((inp_info%bend_a(k))==(name_atm_a)) .AND. & 01097 ((inp_info%bend_b(k))==(name_atm_b)) .AND. & 01098 ((inp_info%bend_c(k))==(name_atm_c))) .OR. & 01099 (((inp_info%bend_a(k))==(name_atm_c)) .AND. & 01100 ((inp_info%bend_b(k))==(name_atm_b)) .AND. & 01101 ((inp_info%bend_c(k))==(name_atm_a))) ) THEN 01102 bend_list(j)%bend_kind%id_type = inp_info%bend_kind(k) 01103 bend_list(j)%bend_kind%k = inp_info%bend_k(k) 01104 bend_list(j)%bend_kind%theta0 = inp_info%bend_theta0(k) 01105 bend_list(j)%bend_kind%cb = inp_info%bend_cb(k) 01106 bend_list(j)%bend_kind%r012 = inp_info%bend_r012(k) 01107 bend_list(j)%bend_kind%r032 = inp_info%bend_r032(k) 01108 bend_list(j)%bend_kind%kbs12 = inp_info%bend_kbs12(k) 01109 bend_list(j)%bend_kind%kbs32 = inp_info%bend_kbs32(k) 01110 bend_list(j)%bend_kind%kss = inp_info%bend_kss(k) 01111 CALL issue_duplications(found,routineP, "Bend", name_atm_a, name_atm_b,& 01112 name_atm_c, error=error) 01113 found = .TRUE. 01114 EXIT 01115 END IF 01116 END DO 01117 END IF 01118 01119 IF(.NOT.found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 01120 atm2=TRIM(name_atm_b),& 01121 atm3=TRIM(name_atm_c),& 01122 fatal=fatal,& 01123 type_name="Angle",& 01124 array=Ainfo,& 01125 error=error) 01126 ! QM/MM modifications 01127 IF (only_qm) THEN 01128 bend_list(j)%id_type = do_ff_undef 01129 bend_list(j)%bend_kind%id_type = do_ff_undef 01130 END IF 01131 END DO 01132 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01133 bend_list=bend_list) 01134 END DO 01135 CALL timestop(handle2) 01136 01137 END SUBROUTINE force_field_pack_bend 01138 01139 ! ***************************************************************************** 01142 SUBROUTINE force_field_pack_ub (particle_set, molecule_kind_set, molecule_set, & 01143 Ainfo, chm_info, inp_info, iw, error) 01144 01145 TYPE(particle_type), DIMENSION(:), 01146 POINTER :: particle_set 01147 TYPE(molecule_kind_type), DIMENSION(:), 01148 POINTER :: molecule_kind_set 01149 TYPE(molecule_type), DIMENSION(:), 01150 POINTER :: molecule_set 01151 CHARACTER(LEN=default_string_length), 01152 DIMENSION(:), POINTER :: Ainfo 01153 TYPE(charmm_info_type), POINTER :: chm_info 01154 TYPE(input_info_type), POINTER :: inp_info 01155 INTEGER :: iw 01156 TYPE(cp_error_type), INTENT(INOUT) :: error 01157 01158 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_ub', 01159 routineP = moduleN//':'//routineN 01160 01161 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, 01162 name_atm_c 01163 INTEGER :: atm_a, atm_b, atm_c, first, 01164 handle2, i, j, k, last, 01165 natom, nub 01166 INTEGER, DIMENSION(:), POINTER :: molecule_list 01167 LOGICAL :: failure, found, only_qm 01168 TYPE(atomic_kind_type), POINTER :: atomic_kind 01169 TYPE(molecule_kind_type), POINTER :: molecule_kind 01170 TYPE(molecule_type), POINTER :: molecule 01171 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 01172 01173 failure = .FALSE. 01174 01175 CALL timeset(routineN,handle2) 01176 DO i=1,SIZE(molecule_kind_set) 01177 molecule_kind => molecule_kind_set(i) 01178 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01179 molecule_list=molecule_list,& 01180 natom=natom,& 01181 nub=nub,ub_list=ub_list) 01182 molecule=>molecule_set(molecule_list(1)) 01183 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 01184 DO j=1,nub 01185 atm_a = ub_list(j)%a 01186 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 01187 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01188 name=name_atm_a) 01189 atm_b = ub_list(j)%b 01190 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 01191 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01192 name=name_atm_b) 01193 atm_c = ub_list(j)%c 01194 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 01195 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01196 name=name_atm_c) 01197 found = .FALSE. 01198 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c) 01199 CALL uppercase(name_atm_a) 01200 CALL uppercase(name_atm_b) 01201 CALL uppercase(name_atm_c) 01202 01203 ! loop over params from GROMOS 01204 ! ikuo - None that I know... 01205 01206 ! loop over params from CHARMM 01207 IF(ASSOCIATED(chm_info%ub_a)) THEN 01208 DO k=1,SIZE(chm_info%ub_a) 01209 IF( (((chm_info%ub_a(k))==(name_atm_a)) .AND. & 01210 ((chm_info%ub_b(k))==(name_atm_b)) .AND. & 01211 ((chm_info%ub_c(k))==(name_atm_c))) .OR. & 01212 (((chm_info%ub_a(k))==(name_atm_c)) .AND. & 01213 ((chm_info%ub_b(k))==(name_atm_b)) .AND. & 01214 ((chm_info%ub_c(k))==(name_atm_a))) ) THEN 01215 ub_list(j)%ub_kind%id_type = do_ff_charmm 01216 ub_list(j)%ub_kind%k(1) = chm_info%ub_k(k) 01217 ub_list(j)%ub_kind%r0 = chm_info%ub_r0(k) 01218 IF(iw>0) WRITE(iw,*) " Found UB ",TRIM(name_atm_a)," ",& 01219 TRIM(name_atm_b)," ",TRIM(name_atm_c) 01220 CALL issue_duplications(found,routineP, "Urey-Bradley", name_atm_a,& 01221 name_atm_b, name_atm_c, error=error) 01222 found = .TRUE. 01223 EXIT 01224 END IF 01225 END DO 01226 END IF 01227 01228 ! loop over params from AMBER 01229 ! teo - None that I know... 01230 01231 ! always have the input param last to overwrite all the other ones 01232 IF(ASSOCIATED(inp_info%ub_a)) THEN 01233 DO k=1,SIZE(inp_info%ub_a) 01234 IF( (((inp_info%ub_a(k))==(name_atm_a)) .AND. & 01235 ((inp_info%ub_b(k))==(name_atm_b)) .AND. & 01236 ((inp_info%ub_c(k))==(name_atm_c))) .OR. & 01237 (((inp_info%ub_a(k))==(name_atm_c)) .AND. & 01238 ((inp_info%ub_b(k))==(name_atm_b)) .AND. & 01239 ((inp_info%ub_c(k))==(name_atm_a))) ) THEN 01240 ub_list(j)%ub_kind%id_type = inp_info%ub_kind(k) 01241 ub_list(j)%ub_kind%k(:) = inp_info%ub_k(:,k) 01242 ub_list(j)%ub_kind%r0 = inp_info%ub_r0(k) 01243 CALL issue_duplications(found,routineP, "Urey-Bradley", name_atm_a,& 01244 name_atm_b, name_atm_c, error=error) 01245 found = .TRUE. 01246 EXIT 01247 END IF 01248 END DO 01249 END IF 01250 01251 IF(.NOT.found) THEN 01252 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 01253 atm2=TRIM(name_atm_b),& 01254 atm3=TRIM(name_atm_c),& 01255 type_name="Urey-Bradley",& 01256 array=Ainfo,& 01257 error=error) 01258 ub_list(j)%id_type = do_ff_undef 01259 ub_list(j)%ub_kind%id_type = do_ff_undef 01260 ub_list(j)%ub_kind%k = 0.0_dp 01261 ub_list(j)%ub_kind%r0 = 0.0_dp 01262 END IF 01263 01264 ! QM/MM modifications 01265 IF (only_qm) THEN 01266 ub_list(j)%id_type = do_ff_undef 01267 ub_list(j)%ub_kind%id_type = do_ff_undef 01268 END IF 01269 END DO 01270 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01271 ub_list=ub_list) 01272 END DO 01273 CALL timestop(handle2) 01274 01275 END SUBROUTINE force_field_pack_ub 01276 01277 ! ***************************************************************************** 01280 SUBROUTINE force_field_pack_tors (particle_set, molecule_kind_set, molecule_set, & 01281 Ainfo, chm_info, inp_info, gro_info, amb_info, iw, error) 01282 01283 TYPE(particle_type), DIMENSION(:), 01284 POINTER :: particle_set 01285 TYPE(molecule_kind_type), DIMENSION(:), 01286 POINTER :: molecule_kind_set 01287 TYPE(molecule_type), DIMENSION(:), 01288 POINTER :: molecule_set 01289 CHARACTER(LEN=default_string_length), 01290 DIMENSION(:), POINTER :: Ainfo 01291 TYPE(charmm_info_type), POINTER :: chm_info 01292 TYPE(input_info_type), POINTER :: inp_info 01293 TYPE(gromos_info_type), POINTER :: gro_info 01294 TYPE(amber_info_type), POINTER :: amb_info 01295 INTEGER, INTENT(IN) :: iw 01296 TYPE(cp_error_type), INTENT(INOUT) :: error 01297 01298 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_tors', 01299 routineP = moduleN//':'//routineN 01300 01301 CHARACTER(LEN=default_string_length) :: ldum, name_atm_a, name_atm_b, 01302 name_atm_c, name_atm_d 01303 INTEGER :: atm_a, atm_b, atm_c, atm_d, 01304 first, handle2, i, imul, 01305 itype, j, k, last, natom, 01306 ntorsion 01307 INTEGER, DIMENSION(:), POINTER :: molecule_list 01308 LOGICAL :: failure, found, only_qm 01309 TYPE(atomic_kind_type), POINTER :: atomic_kind 01310 TYPE(molecule_kind_type), POINTER :: molecule_kind 01311 TYPE(molecule_type), POINTER :: molecule 01312 TYPE(torsion_type), DIMENSION(:), 01313 POINTER :: torsion_list 01314 01315 CALL timeset(routineN,handle2) 01316 failure = .FALSE. 01317 01318 DO i=1,SIZE(molecule_kind_set) 01319 molecule_kind => molecule_kind_set(i) 01320 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01321 molecule_list=molecule_list,& 01322 natom=natom,& 01323 ntorsion=ntorsion,torsion_list=torsion_list) 01324 molecule=>molecule_set(molecule_list(1)) 01325 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 01326 DO j=1,ntorsion 01327 IF(torsion_list(j)%torsion_kind%id_type == do_ff_undef) THEN 01328 atm_a = torsion_list(j)%a 01329 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 01330 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01331 name=name_atm_a) 01332 atm_b = torsion_list(j)%b 01333 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 01334 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01335 name=name_atm_b) 01336 atm_c = torsion_list(j)%c 01337 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 01338 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01339 name=name_atm_c) 01340 atm_d = torsion_list(j)%d 01341 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 01342 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01343 name=name_atm_d) 01344 found = .FALSE. 01345 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 01346 CALL uppercase(name_atm_a) 01347 CALL uppercase(name_atm_b) 01348 CALL uppercase(name_atm_c) 01349 CALL uppercase(name_atm_d) 01350 01351 ! loop over params from GROMOS 01352 IF(ASSOCIATED(gro_info%torsion_k)) THEN 01353 k=SIZE(gro_info%torsion_k) 01354 itype = torsion_list(j)%itype 01355 IF(itype>0) THEN 01356 CALL reallocate(torsion_list(j)%torsion_kind%k,1,1) 01357 CALL reallocate(torsion_list(j)%torsion_kind%m,1,1) 01358 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1) 01359 torsion_list(j)%torsion_kind%nmul = 1 01360 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype) 01361 torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype) 01362 torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype) 01363 ELSE 01364 CALL reallocate(torsion_list(j)%torsion_kind%k,1,1) 01365 CALL reallocate(torsion_list(j)%torsion_kind%m,1,1) 01366 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,1) 01367 torsion_list(j)%torsion_kind%nmul = 1 01368 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype/k) 01369 torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype/k) 01370 torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype/k) 01371 END IF 01372 torsion_list(j)%torsion_kind%id_type = gro_info%ff_gromos_type 01373 torsion_list(j)%id_type = gro_info%ff_gromos_type 01374 found = .TRUE. 01375 imul = torsion_list(j)%torsion_kind%nmul 01376 END IF 01377 01378 ! loop over params from CHARMM 01379 IF(ASSOCIATED(chm_info%torsion_a)) THEN 01380 DO k=1,SIZE(chm_info%torsion_a) 01381 IF( (((chm_info%torsion_a(k))==(name_atm_a)) .AND. & 01382 ((chm_info%torsion_b(k))==(name_atm_b)) .AND. & 01383 ((chm_info%torsion_c(k))==(name_atm_c)) .AND. & 01384 ((chm_info%torsion_d(k))==(name_atm_d))) .OR. & 01385 (((chm_info%torsion_a(k))==(name_atm_d)) .AND. & 01386 ((chm_info%torsion_b(k))==(name_atm_c)) .AND. & 01387 ((chm_info%torsion_c(k))==(name_atm_b)) .AND. & 01388 ((chm_info%torsion_d(k))==(name_atm_a))) ) THEN 01389 imul = torsion_list(j)%torsion_kind%nmul + 1 01390 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) 01391 CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) 01392 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) 01393 torsion_list(j)%torsion_kind%id_type = do_ff_charmm 01394 torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) 01395 torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) 01396 torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) 01397 torsion_list(j)%torsion_kind%nmul = imul 01398 found = .TRUE. 01399 END IF 01400 END DO 01401 01402 IF(.NOT.found) THEN 01403 DO k=1,SIZE(chm_info%torsion_a) 01404 IF( (((chm_info%torsion_a(k))==("X")) .AND. & 01405 ((chm_info%torsion_b(k))==(name_atm_b)) .AND. & 01406 ((chm_info%torsion_c(k))==(name_atm_c)) .AND. & 01407 ((chm_info%torsion_d(k))==("X"))) .OR. & 01408 (((chm_info%torsion_a(k))==("X")) .AND. & 01409 ((chm_info%torsion_b(k))==(name_atm_c)) .AND. & 01410 ((chm_info%torsion_c(k))==(name_atm_b)) .AND. & 01411 ((chm_info%torsion_d(k))==("X"))) ) THEN 01412 imul = torsion_list(j)%torsion_kind%nmul + 1 01413 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) 01414 CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) 01415 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) 01416 torsion_list(j)%torsion_kind%id_type = do_ff_charmm 01417 torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) 01418 torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) 01419 torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) 01420 torsion_list(j)%torsion_kind%nmul = imul 01421 found = .TRUE. 01422 END IF 01423 END DO 01424 END IF 01425 END IF 01426 01427 ! loop over params from AMBER 01428 IF(ASSOCIATED(amb_info%torsion_a)) THEN 01429 DO k=1,SIZE(amb_info%torsion_a) 01430 IF( (((amb_info%torsion_a(k))==(name_atm_a)) .AND. & 01431 ((amb_info%torsion_b(k))==(name_atm_b)) .AND. & 01432 ((amb_info%torsion_c(k))==(name_atm_c)) .AND. & 01433 ((amb_info%torsion_d(k))==(name_atm_d))) .OR. & 01434 (((amb_info%torsion_a(k))==(name_atm_d)) .AND. & 01435 ((amb_info%torsion_b(k))==(name_atm_c)) .AND. & 01436 ((amb_info%torsion_c(k))==(name_atm_b)) .AND. & 01437 ((amb_info%torsion_d(k))==(name_atm_a))) ) THEN 01438 imul = torsion_list(j)%torsion_kind%nmul + 1 01439 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) 01440 CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) 01441 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) 01442 torsion_list(j)%torsion_kind%id_type = do_ff_amber 01443 torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) 01444 torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) 01445 torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) 01446 torsion_list(j)%torsion_kind%nmul = imul 01447 found = .TRUE. 01448 END IF 01449 END DO 01450 01451 IF(.NOT.found) THEN 01452 DO k=1,SIZE(amb_info%torsion_a) 01453 IF( (((amb_info%torsion_a(k))==("X")) .AND. & 01454 ((amb_info%torsion_b(k))==(name_atm_b)) .AND. & 01455 ((amb_info%torsion_c(k))==(name_atm_c)) .AND. & 01456 ((amb_info%torsion_d(k))==("X"))) .OR. & 01457 (((amb_info%torsion_a(k))==("X")) .AND. & 01458 ((amb_info%torsion_b(k))==(name_atm_c)) .AND. & 01459 ((amb_info%torsion_c(k))==(name_atm_b)) .AND. & 01460 ((amb_info%torsion_d(k))==("X"))) ) THEN 01461 imul = torsion_list(j)%torsion_kind%nmul + 1 01462 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) 01463 CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) 01464 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) 01465 torsion_list(j)%torsion_kind%id_type = do_ff_amber 01466 torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) 01467 torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) 01468 torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) 01469 torsion_list(j)%torsion_kind%nmul = imul 01470 found = .TRUE. 01471 END IF 01472 END DO 01473 END IF 01474 END IF 01475 01476 ! always have the input param last to overwrite all the other ones 01477 IF(ASSOCIATED(inp_info%torsion_a)) THEN 01478 DO k=1,SIZE(inp_info%torsion_a) 01479 IF( (((inp_info%torsion_a(k))==(name_atm_a)) .AND. & 01480 ((inp_info%torsion_b(k))==(name_atm_b)) .AND. & 01481 ((inp_info%torsion_c(k))==(name_atm_c)) .AND. & 01482 ((inp_info%torsion_d(k))==(name_atm_d))) .OR. & 01483 (((inp_info%torsion_a(k))==(name_atm_d)) .AND. & 01484 ((inp_info%torsion_b(k))==(name_atm_c)) .AND. & 01485 ((inp_info%torsion_c(k))==(name_atm_b)) .AND. & 01486 ((inp_info%torsion_d(k))==(name_atm_a))) ) THEN 01487 imul = torsion_list(j)%torsion_kind%nmul + 1 01488 CALL reallocate(torsion_list(j)%torsion_kind%k,1,imul) 01489 CALL reallocate(torsion_list(j)%torsion_kind%m,1,imul) 01490 CALL reallocate(torsion_list(j)%torsion_kind%phi0,1,imul) 01491 torsion_list(j)%torsion_kind%id_type = inp_info%torsion_kind(k) 01492 torsion_list(j)%torsion_kind%k(imul) = inp_info%torsion_k(k) 01493 torsion_list(j)%torsion_kind%m(imul) = inp_info%torsion_m(k) 01494 torsion_list(j)%torsion_kind%phi0(imul) = inp_info%torsion_phi0(k) 01495 torsion_list(j)%torsion_kind%nmul = imul 01496 found = .TRUE. 01497 END IF 01498 END DO 01499 END IF 01500 01501 IF(.NOT.found) THEN 01502 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 01503 atm2=TRIM(name_atm_b),& 01504 atm3=TRIM(name_atm_c),& 01505 atm4=TRIM(name_atm_d),& 01506 type_name="Torsion",& 01507 array=Ainfo,& 01508 error=error) 01509 torsion_list(j)%torsion_kind%id_type = do_ff_undef 01510 torsion_list(j)%id_type = do_ff_undef 01511 ELSE 01512 ldum = cp_to_string(imul) 01513 IF ((imul /= 1).AND.(iw>0))& 01514 WRITE(iw,'(/,2("UTIL_INFO| ",A,/))')& 01515 "Multiple Torsion declarations: "//TRIM(name_atm_a)//& 01516 ","//TRIM(name_atm_b)//","//TRIM(name_atm_c)//" and "//TRIM(name_atm_d),& 01517 "Number of defined torsions: "//TRIM(ldum)//" ." 01518 END IF 01519 ! 01520 ! QM/MM modifications 01521 ! 01522 IF (only_qm) THEN 01523 IF (iw>0) WRITE(iw,*)" Torsion PARAM between QM atoms ",j," : ",& 01524 TRIM(name_atm_a)," ",& 01525 TRIM(name_atm_b)," ",& 01526 TRIM(name_atm_c)," ",& 01527 TRIM(name_atm_d)," ",& 01528 torsion_list(j)%a,& 01529 torsion_list(j)%b,& 01530 torsion_list(j)%c,& 01531 torsion_list(j)%d 01532 torsion_list(j)%torsion_kind%id_type = do_ff_undef 01533 torsion_list(j)%id_type = do_ff_undef 01534 END IF 01535 END IF 01536 END DO 01537 CALL set_molecule_kind(molecule_kind=molecule_kind,& 01538 torsion_list=torsion_list) 01539 END DO 01540 CALL timestop(handle2) 01541 01542 END SUBROUTINE force_field_pack_tors 01543 01544 ! ***************************************************************************** 01547 SUBROUTINE force_field_pack_impr (particle_set, molecule_kind_set, molecule_set, & 01548 Ainfo, chm_info, inp_info, gro_info, amb_info, error) 01549 01550 TYPE(particle_type), DIMENSION(:), 01551 POINTER :: particle_set 01552 TYPE(molecule_kind_type), DIMENSION(:), 01553 POINTER :: molecule_kind_set 01554 TYPE(molecule_type), DIMENSION(:), 01555 POINTER :: molecule_set 01556 CHARACTER(LEN=default_string_length), 01557 DIMENSION(:), POINTER :: Ainfo 01558 TYPE(charmm_info_type), POINTER :: chm_info 01559 TYPE(input_info_type), POINTER :: inp_info 01560 TYPE(gromos_info_type), POINTER :: gro_info 01561 TYPE(amber_info_type), POINTER :: amb_info 01562 TYPE(cp_error_type), INTENT(INOUT) :: error 01563 01564 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_impr', 01565 routineP = moduleN//':'//routineN 01566 01567 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, 01568 name_atm_c, name_atm_d 01569 INTEGER :: atm_a, atm_b, atm_c, atm_d, 01570 first, handle2, i, itype, j, 01571 k, last, natom, nimpr 01572 INTEGER, DIMENSION(:), POINTER :: molecule_list 01573 LOGICAL :: failure, found, only_qm 01574 TYPE(atomic_kind_type), POINTER :: atomic_kind 01575 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 01576 TYPE(molecule_kind_type), POINTER :: molecule_kind 01577 TYPE(molecule_type), POINTER :: molecule 01578 01579 CALL timeset(routineN,handle2) 01580 failure = .FALSE. 01581 01582 DO i=1,SIZE(molecule_kind_set) 01583 molecule_kind => molecule_kind_set(i) 01584 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01585 molecule_list=molecule_list,& 01586 natom=natom,& 01587 nimpr=nimpr,impr_list=impr_list) 01588 molecule=>molecule_set(molecule_list(1)) 01589 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 01590 DO j=1,nimpr 01591 atm_a = impr_list(j)%a 01592 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 01593 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01594 name=name_atm_a) 01595 atm_b = impr_list(j)%b 01596 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 01597 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01598 name=name_atm_b) 01599 atm_c = impr_list(j)%c 01600 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 01601 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01602 name=name_atm_c) 01603 atm_d = impr_list(j)%d 01604 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 01605 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01606 name=name_atm_d) 01607 found = .FALSE. 01608 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 01609 CALL uppercase(name_atm_a) 01610 CALL uppercase(name_atm_b) 01611 CALL uppercase(name_atm_c) 01612 CALL uppercase(name_atm_d) 01613 01614 ! loop over params from GROMOS 01615 IF(ASSOCIATED(gro_info%impr_k)) THEN 01616 k=SIZE(gro_info%impr_k) 01617 itype = impr_list(j)%itype 01618 IF(itype>0) THEN 01619 impr_list(j)%impr_kind%k = gro_info%impr_k(itype) 01620 impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype) 01621 ELSE 01622 impr_list(j)%impr_kind%k = gro_info%impr_k(itype) 01623 impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype) 01624 END IF 01625 found = .TRUE. 01626 impr_list(j)%impr_kind%id_type = gro_info%ff_gromos_type 01627 impr_list(j)%id_type = gro_info%ff_gromos_type 01628 END IF 01629 01630 ! loop over params from CHARMM 01631 IF(ASSOCIATED(chm_info%impr_a)) THEN 01632 DO k=1,SIZE(chm_info%impr_a) 01633 IF( (((chm_info%impr_a(k))==(name_atm_a)) .AND. & 01634 ((chm_info%impr_b(k))==(name_atm_b)) .AND. & 01635 ((chm_info%impr_c(k))==(name_atm_c)) .AND. & 01636 ((chm_info%impr_d(k))==(name_atm_d))) .OR. & 01637 (((chm_info%impr_a(k))==(name_atm_d)) .AND. & 01638 ((chm_info%impr_b(k))==(name_atm_c)) .AND. & 01639 ((chm_info%impr_c(k))==(name_atm_b)) .AND. & 01640 ((chm_info%impr_d(k))==(name_atm_a))) ) THEN 01641 impr_list(j)%impr_kind%id_type = do_ff_charmm 01642 impr_list(j)%impr_kind%k = chm_info%impr_k(k) 01643 impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k) 01644 CALL issue_duplications(found,routineP, "Impropers", name_atm_a, name_atm_b,& 01645 name_atm_c, name_atm_d, error=error) 01646 found = .TRUE. 01647 EXIT 01648 END IF 01649 END DO 01650 IF(.NOT.found) THEN 01651 DO k=1,SIZE(chm_info%impr_a) 01652 IF( (((chm_info%impr_a(k))==(name_atm_a)) .AND. & 01653 ((chm_info%impr_b(k))==("X")) .AND. & 01654 ((chm_info%impr_c(k))==("X")) .AND. & 01655 ((chm_info%impr_d(k))==(name_atm_d))) .OR. & 01656 (((chm_info%impr_a(k))==(name_atm_d)) .AND. & 01657 ((chm_info%impr_b(k))==("X")) .AND. & 01658 ((chm_info%impr_c(k))==("X")) .AND. & 01659 ((chm_info%impr_d(k))==(name_atm_a))) ) THEN 01660 impr_list(j)%impr_kind%id_type = do_ff_charmm 01661 impr_list(j)%impr_kind%k = chm_info%impr_k(k) 01662 impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k) 01663 CALL issue_duplications(found,routineP, "Impropers", name_atm_a, name_atm_b,& 01664 name_atm_c,name_atm_d, error=error) 01665 found = .TRUE. 01666 EXIT 01667 END IF 01668 END DO 01669 END IF 01670 END IF 01671 01672 ! loop over params from AMBER not needed since impropers in AMBER 01673 ! are treated like standard torsions 01674 01675 ! always have the input param last to overwrite all the other ones 01676 IF(ASSOCIATED(inp_info%impr_a)) THEN 01677 DO k=1,SIZE(inp_info%impr_a) 01678 IF( ((inp_info%impr_a(k))==(name_atm_a)) .AND. & 01679 ((inp_info%impr_b(k))==(name_atm_b)) .AND. & 01680 ((((inp_info%impr_c(k))==(name_atm_c)) .AND. & 01681 ((inp_info%impr_d(k))==(name_atm_d))) .OR. & 01682 (((inp_info%impr_c(k))==(name_atm_d)) .AND. & 01683 ((inp_info%impr_d(k))==(name_atm_c)))) ) THEN 01684 impr_list(j)%impr_kind%id_type = inp_info%impr_kind(k) 01685 impr_list(j)%impr_kind%k = inp_info%impr_k(k) 01686 IF ( ((inp_info%impr_c(k))==(name_atm_c)) .AND. & 01687 ((inp_info%impr_d(k))==(name_atm_d)) ) THEN 01688 impr_list(j)%impr_kind%phi0 = inp_info%impr_phi0(k) 01689 ELSE 01690 impr_list(j)%impr_kind%phi0 = -inp_info%impr_phi0(k) 01691 ! alternative solutions: 01692 ! - swap impr_list(j)%c with impr_list(j)%d and 01693 ! name_atom_c with name_atom_d and 01694 ! atm_c with atm_d 01695 ! - introduce impr_list(j)%impr_kind%sign. if one, the 01696 ! sign of phi is not changed in mol_force.f90. if minus 01697 ! one, the sign of phi is changed in mol_force.f90 01698 ! similar problems with parameters from charmm pot file 01699 ! above 01700 END IF 01701 CALL issue_duplications(found,routineP, "Impropers", name_atm_a, name_atm_b,& 01702 name_atm_c,name_atm_d, error=error) 01703 found = .TRUE. 01704 EXIT 01705 END IF 01706 END DO 01707 END IF 01708 01709 IF(.NOT.found) THEN 01710 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 01711 atm2=TRIM(name_atm_b),& 01712 atm3=TRIM(name_atm_c),& 01713 atm4=TRIM(name_atm_d),& 01714 type_name="Improper",& 01715 array=Ainfo,& 01716 error=error) 01717 impr_list(j)%impr_kind%k = 0.0_dp 01718 impr_list(j)%impr_kind%phi0 = 0.0_dp 01719 impr_list(j)%impr_kind%id_type = do_ff_undef 01720 impr_list(j)%id_type = do_ff_undef 01721 END IF 01722 ! 01723 ! QM/MM modifications 01724 ! 01725 IF (only_qm) THEN 01726 impr_list(j)%impr_kind%id_type = do_ff_undef 01727 impr_list(j)%id_type = do_ff_undef 01728 END IF 01729 01730 END DO 01731 CALL set_molecule_kind(molecule_kind=molecule_kind,impr_list=impr_list) 01732 END DO 01733 CALL timestop(handle2) 01734 01735 END SUBROUTINE force_field_pack_impr 01736 01737 ! ***************************************************************************** 01743 SUBROUTINE force_field_pack_opbend (particle_set, molecule_kind_set, molecule_set, & 01744 Ainfo, inp_info, error) 01745 01746 TYPE(particle_type), DIMENSION(:), 01747 POINTER :: particle_set 01748 TYPE(molecule_kind_type), DIMENSION(:), 01749 POINTER :: molecule_kind_set 01750 TYPE(molecule_type), DIMENSION(:), 01751 POINTER :: molecule_set 01752 CHARACTER(LEN=default_string_length), 01753 DIMENSION(:), POINTER :: Ainfo 01754 TYPE(input_info_type), POINTER :: inp_info 01755 TYPE(cp_error_type), INTENT(INOUT) :: error 01756 01757 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_opbend', 01758 routineP = moduleN//':'//routineN 01759 01760 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, 01761 name_atm_c, name_atm_d 01762 INTEGER :: atm_a, atm_b, atm_c, atm_d, 01763 first, handle2, i, j, k, 01764 last, natom, nopbend 01765 INTEGER, DIMENSION(:), POINTER :: molecule_list 01766 LOGICAL :: failure, found, only_qm 01767 TYPE(atomic_kind_type), POINTER :: atomic_kind 01768 TYPE(molecule_kind_type), POINTER :: molecule_kind 01769 TYPE(molecule_type), POINTER :: molecule 01770 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 01771 01772 CALL timeset(routineN,handle2) 01773 failure = .FALSE. 01774 01775 DO i=1,SIZE(molecule_kind_set) 01776 molecule_kind => molecule_kind_set(i) 01777 CALL get_molecule_kind(molecule_kind=molecule_kind,& 01778 molecule_list=molecule_list,& 01779 natom=natom,& 01780 nopbend=nopbend,opbend_list=opbend_list) 01781 molecule=>molecule_set(molecule_list(1)) 01782 01783 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 01784 DO j=1,nopbend 01785 atm_a = opbend_list(j)%a 01786 atomic_kind => particle_set(atm_a+first-1)%atomic_kind 01787 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01788 name=name_atm_a) 01789 atm_b = opbend_list(j)%b 01790 atomic_kind => particle_set(atm_b+first-1)%atomic_kind 01791 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01792 name=name_atm_b) 01793 atm_c = opbend_list(j)%c 01794 atomic_kind => particle_set(atm_c+first-1)%atomic_kind 01795 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01796 name=name_atm_c) 01797 atm_d = opbend_list(j)%d 01798 atomic_kind => particle_set(atm_d+first-1)%atomic_kind 01799 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01800 name=name_atm_d) 01801 found = .FALSE. 01802 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 01803 CALL uppercase(name_atm_a) 01804 CALL uppercase(name_atm_b) 01805 CALL uppercase(name_atm_c) 01806 CALL uppercase(name_atm_d) 01807 01808 ! always have the input param last to overwrite all the other ones 01809 IF(ASSOCIATED(inp_info%opbend_a)) THEN 01810 DO k=1,SIZE(inp_info%opbend_a) 01811 IF( ((inp_info%opbend_a(k))==(name_atm_a)) .AND. & 01812 ((inp_info%opbend_d(k))==(name_atm_d)) .AND. & 01813 ((((inp_info%opbend_c(k))==(name_atm_c)) .AND. & 01814 ((inp_info%opbend_b(k))==(name_atm_b))) .OR. & 01815 (((inp_info%opbend_c(k))==(name_atm_b)) .AND. & 01816 ((inp_info%opbend_b(k))==(name_atm_c)))) ) THEN 01817 opbend_list(j)%opbend_kind%id_type = inp_info%opbend_kind(k) 01818 opbend_list(j)%opbend_kind%k = inp_info%opbend_k(k) 01819 IF ( ((inp_info%opbend_c(k))==(name_atm_c)) .AND. & 01820 ((inp_info%opbend_b(k))==(name_atm_b)) ) THEN 01821 opbend_list(j)%opbend_kind%phi0 = inp_info%opbend_phi0(k) 01822 ELSE 01823 opbend_list(j)%opbend_kind%phi0 = -inp_info%opbend_phi0(k) 01824 ! alternative solutions: 01825 ! - swap opbend_list(j)%c with opbend_list(j)%b and 01826 ! name_atom_c with name_atom_b and 01827 ! atm_c with atm_b 01828 ! - introduce opbend_list(j)%opbend_kind%sign. if one, the 01829 ! sign of phi is not changed in mol_force.f90. if minus 01830 ! one, the sign of phi is changed in mol_force.f90 01831 01832 END IF 01833 CALL issue_duplications(found,routineP, "Out of plane bend", name_atm_a, name_atm_b,& 01834 name_atm_c,name_atm_d, error=error) 01835 found = .TRUE. 01836 EXIT 01837 END IF 01838 END DO 01839 END IF 01840 01841 IF(.NOT.found) THEN 01842 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 01843 atm2=TRIM(name_atm_b),& 01844 atm3=TRIM(name_atm_c),& 01845 atm4=TRIM(name_atm_d),& 01846 type_name="Out of plane bend",& 01847 array=Ainfo,& 01848 error=error) 01849 opbend_list(j)%opbend_kind%k = 0.0_dp 01850 opbend_list(j)%opbend_kind%phi0 = 0.0_dp 01851 opbend_list(j)%opbend_kind%id_type = do_ff_undef 01852 opbend_list(j)%id_type = do_ff_undef 01853 END IF 01854 ! 01855 ! QM/MM modifications 01856 ! 01857 IF (only_qm) THEN 01858 opbend_list(j)%opbend_kind%id_type = do_ff_undef 01859 opbend_list(j)%id_type = do_ff_undef 01860 END IF 01861 01862 END DO 01863 CALL set_molecule_kind(molecule_kind=molecule_kind,opbend_list=opbend_list) 01864 END DO 01865 CALL timestop(handle2) 01866 01867 END SUBROUTINE force_field_pack_opbend 01868 01869 ! ***************************************************************************** 01874 SUBROUTINE force_field_pack_charges(charges, charges_section, particle_set, & 01875 my_qmmm, qmmm_env, inp_info, iw4, error) 01876 01877 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 01878 TYPE(section_vals_type), POINTER :: charges_section 01879 TYPE(particle_type), DIMENSION(:), 01880 POINTER :: particle_set 01881 LOGICAL :: my_qmmm 01882 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 01883 TYPE(input_info_type), POINTER :: inp_info 01884 INTEGER :: iw4 01885 TYPE(cp_error_type), INTENT(INOUT) :: error 01886 01887 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charges', 01888 routineP = moduleN//':'//routineN 01889 01890 CHARACTER(LEN=default_string_length) :: atmname 01891 INTEGER :: handle, iatom, ilink, j, 01892 nval, stat 01893 LOGICAL :: failure, found_p, 01894 is_link_atom, is_ok, 01895 only_manybody, only_qm 01896 REAL(KIND=dp) :: charge, charge_tot, rval, 01897 scale_factor 01898 TYPE(atomic_kind_type), POINTER :: atomic_kind 01899 TYPE(cp_sll_val_type), POINTER :: list 01900 TYPE(fist_potential_type), POINTER :: fist_potential 01901 TYPE(val_type), POINTER :: val 01902 01903 failure = .FALSE. 01904 01905 CALL timeset(routineN,handle) 01906 charge_tot = 0.0_dp 01907 NULLIFY(list) 01908 01909 ! Not implemented for core-shell 01910 IF(ASSOCIATED(inp_info%shell_list)) THEN 01911 CALL cp_unimplemented_error(fromWhere=routineP, & 01912 message="Array of charges not implemented for CORE-SHELL model!!",& 01913 error=error, error_level=cp_failure_level) 01914 END IF 01915 01916 ! Allocate array to particle_set size 01917 CPPostcondition(.NOT.(ASSOCIATED(charges)),cp_failure_level,routineP,error,failure) 01918 ALLOCATE(charges(SIZE(particle_set)),stat=stat) 01919 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01920 01921 ! Fill with input values 01922 CALL section_vals_val_get(charges_section,"_DEFAULT_KEYWORD_",n_rep_val=nval,error=error) 01923 CPPostcondition(nval==SIZE(charges),cp_failure_level,routineP,error,failure) 01924 CALL section_vals_list_get(charges_section,"_DEFAULT_KEYWORD_",list=list,error=error) 01925 DO iatom=1,nval 01926 ! we use only the first default_string_length characters of each line 01927 is_ok=cp_sll_val_next(list,val,error=error) 01928 CALL val_get(val,r_val=rval,error=error) 01929 ! assign values 01930 charges(iatom)=rval 01931 01932 ! Perform a post-processing 01933 atomic_kind => particle_set(iatom)%atomic_kind 01934 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01935 fist_potential=fist_potential,& 01936 name=atmname) 01937 CALL get_potential(potential=fist_potential, qeff=charge) 01938 01939 only_qm = qmmm_ff_precond_only_qm(id1=atmname,is_link=is_link_atom) 01940 CALL uppercase(atmname) 01941 CALL cp_assert(charge==-HUGE(0.0_dp),cp_warning_level,cp_assertion_failed,routineP,& 01942 "The charge for atom index ("//cp_to_string(iatom)//") and atom name ("//& 01943 TRIM(atmname)//") was already defined. The charge associated to this kind"//& 01944 " will be set to an uninitialized value and only the atom specific charge will be used! "//& 01945 CPSourceFileRef,& 01946 only_ionode=.TRUE.) 01947 charge=-HUGE(0.0_dp) 01948 01949 ! Check if the potential really requires the charge definition.. 01950 IF (ASSOCIATED(inp_info%nonbonded)) THEN 01951 IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN 01952 ! Let's find the nonbonded potential where this atom is involved 01953 only_manybody = .TRUE. 01954 found_p = .FALSE. 01955 DO j = 1, SIZE(inp_info%nonbonded%pot) 01956 IF ( atmname==inp_info%nonbonded%pot(j)%pot %at1.OR.& 01957 atmname==inp_info%nonbonded%pot(j)%pot %at2) THEN 01958 SELECT CASE(inp_info%nonbonded%pot(j)%pot%type(1)) 01959 CASE (ea_type,tersoff_type,siepmann_type) 01960 ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential 01961 ! Do nothing.. 01962 CASE DEFAULT 01963 only_manybody = .FALSE. 01964 EXIT 01965 END SELECT 01966 found_p = .TRUE. 01967 END IF 01968 END DO 01969 IF (only_manybody.AND.found_p) THEN 01970 charges(iatom)= 0.0_dp 01971 END IF 01972 END IF 01973 END IF 01974 ! 01975 ! QM/MM modifications 01976 ! 01977 IF (only_qm.AND.my_qmmm) THEN 01978 IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN 01979 scale_factor = 0.0_dp 01980 IF (is_link_atom) THEN 01981 ! 01982 ! Find the scaling factor... 01983 ! 01984 DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) 01985 IF (iatom == qmmm_env%mm_link_atoms(ilink)) EXIT 01986 END DO 01987 CPPostcondition(ilink <= SIZE(qmmm_env%mm_link_atoms),cp_failure_level,routineP,error,failure) 01988 scale_factor = qmmm_env%fist_scale_charge_link(ilink) 01989 END IF 01990 charges(iatom) = charges(iatom) * scale_factor 01991 END IF 01992 END IF 01993 END DO 01994 ! Sum up total charges for IO 01995 charge_tot = SUM(charges) 01996 ! Print Total Charge of the system 01997 IF (iw4>0) THEN 01998 WRITE(iw4,FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)')charge_tot 01999 END IF 02000 CALL timestop(handle) 02001 END SUBROUTINE force_field_pack_charges 02002 02003 ! ***************************************************************************** 02007 SUBROUTINE force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4,& 02008 Ainfo, my_qmmm, inp_info, error) 02009 02010 TYPE(atomic_kind_type), DIMENSION(:), 02011 POINTER :: atomic_kind_set 02012 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 02013 LOGICAL :: fatal 02014 INTEGER :: iw, iw4 02015 CHARACTER(LEN=default_string_length), 02016 DIMENSION(:), POINTER :: Ainfo 02017 LOGICAL :: my_qmmm 02018 TYPE(input_info_type), POINTER :: inp_info 02019 TYPE(cp_error_type), INTENT(INOUT) :: error 02020 02021 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charge', 02022 routineP = moduleN//':'//routineN 02023 02024 CHARACTER(LEN=default_string_length) :: atmname 02025 INTEGER :: handle, i, ilink, j 02026 INTEGER, DIMENSION(:), POINTER :: my_atom_list 02027 LOGICAL :: failure, found, found_p, 02028 is_link_atom, is_shell, 02029 only_manybody, only_qm 02030 REAL(KIND=dp) :: charge, charge_tot, 02031 cs_charge, scale_factor 02032 TYPE(atomic_kind_type), POINTER :: atomic_kind 02033 TYPE(fist_potential_type), POINTER :: fist_potential 02034 02035 failure = .FALSE. 02036 02037 CALL timeset(routineN,handle) 02038 charge_tot = 0.0_dp 02039 DO i=1,SIZE(atomic_kind_set) 02040 atomic_kind => atomic_kind_set(i) 02041 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02042 fist_potential=fist_potential,& 02043 atom_list=my_atom_list,& 02044 name=atmname) 02045 CALL get_potential(potential=fist_potential, qeff=charge) 02046 02047 is_shell = .FALSE. 02048 found = .FALSE. 02049 only_qm = qmmm_ff_precond_only_qm(id1=atmname,is_link=is_link_atom) 02050 CALL uppercase(atmname) 02051 IF(charge/=-HUGE(0.0_dp)) found = .TRUE. 02052 02053 ! Always have the input param last to overwrite all the other ones 02054 IF(ASSOCIATED(inp_info%charge_atm)) THEN 02055 DO j=1,SIZE(inp_info%charge_atm) 02056 IF (iw>0) WRITE(iw,*)"Charge Checking ::",TRIM(inp_info%charge_atm(j)),atmname 02057 IF((inp_info%charge_atm(j))==atmname) THEN 02058 charge = inp_info%charge(j) 02059 CALL issue_duplications(found,routineP, "Charge", atmname, error=error) 02060 found = .TRUE. 02061 END IF 02062 END DO 02063 END IF 02064 ! Check if the ATOM type has a core-shell associated.. In this case 02065 ! print a warning: the CHARGE will not be used if defined.. or we can 02066 ! even skip its definition.. 02067 IF(ASSOCIATED(inp_info%shell_list)) THEN 02068 DO j=1,SIZE(inp_info%shell_list) 02069 IF((inp_info%shell_list(j)%atm_name)==atmname) THEN 02070 is_shell = .TRUE. 02071 cs_charge = inp_info%shell_list(j)%shell%charge_core+& 02072 inp_info%shell_list(j)%shell%charge_shell 02073 charge = 0.0_dp 02074 IF (found) THEN 02075 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 02076 "CORE-SHELL model defined for KIND ("//TRIM(atmname)//")"//& 02077 " ignoring charge definition! "//& 02078 CPSourceFileRef,& 02079 only_ionode=.TRUE.) 02080 ELSE 02081 found = .TRUE. 02082 END IF 02083 END IF 02084 END DO 02085 END IF 02086 ! Check if the potential really requires the charge definition.. 02087 IF (ASSOCIATED(inp_info%nonbonded)) THEN 02088 IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN 02089 ! Let's find the nonbonded potential where this atom is involved 02090 only_manybody = .TRUE. 02091 found_p = .FALSE. 02092 DO j = 1, SIZE(inp_info%nonbonded%pot) 02093 IF ( atmname==inp_info%nonbonded%pot(j)%pot %at1.OR.& 02094 atmname==inp_info%nonbonded%pot(j)%pot %at2) THEN 02095 SELECT CASE(inp_info%nonbonded%pot(j)%pot%type(1)) 02096 CASE (ea_type,tersoff_type,siepmann_type) 02097 ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential 02098 ! Do nothing.. 02099 CASE DEFAULT 02100 only_manybody = .FALSE. 02101 EXIT 02102 END SELECT 02103 found_p = .TRUE. 02104 END IF 02105 END DO 02106 IF (only_manybody.AND.found_p) THEN 02107 charge = 0.0_dp 02108 found = .TRUE. 02109 END IF 02110 END IF 02111 END IF 02112 IF (.NOT.found) THEN 02113 ! Set the charge to zero anyway in case the user decides to ignore 02114 ! missing critical parameters. 02115 charge = 0.0_dp 02116 CALL store_FF_missing_par(atm1=TRIM(atmname),& 02117 fatal=fatal,& 02118 type_name="Charge",& 02119 array=Ainfo,& 02120 error=error) 02121 END IF 02122 ! 02123 ! QM/MM modifications 02124 ! 02125 IF (only_qm.AND.my_qmmm) THEN 02126 IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN 02127 scale_factor = 0.0_dp 02128 IF (is_link_atom) THEN 02129 ! 02130 ! Find the scaling factor... 02131 ! 02132 DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) 02133 IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT 02134 END DO 02135 CPPostcondition(ilink <= SIZE(qmmm_env%mm_link_atoms),cp_failure_level,routineP,error,failure) 02136 scale_factor = qmmm_env%fist_scale_charge_link(ilink) 02137 END IF 02138 charge = charge * scale_factor 02139 END IF 02140 END IF 02141 02142 CALL set_potential(potential=fist_potential, qeff=charge) 02143 ! Sum up total charges for IO 02144 IF (found) THEN 02145 IF (is_shell) THEN 02146 charge_tot = charge_tot + atomic_kind%natom*cs_charge 02147 ELSE 02148 charge_tot = charge_tot + atomic_kind%natom*charge 02149 END IF 02150 END IF 02151 END DO 02152 ! Print Total Charge of the system 02153 IF (iw4>0) THEN 02154 WRITE(iw4,FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)')charge_tot 02155 END IF 02156 CALL timestop(handle) 02157 END SUBROUTINE force_field_pack_charge 02158 02159 ! ***************************************************************************** 02163 SUBROUTINE force_field_pack_radius(atomic_kind_set, iw, subsys_section, error) 02164 02165 TYPE(atomic_kind_type), DIMENSION(:), 02166 POINTER :: atomic_kind_set 02167 INTEGER, INTENT(IN) :: iw 02168 TYPE(section_vals_type), POINTER :: subsys_section 02169 TYPE(cp_error_type), INTENT(INOUT) :: error 02170 02171 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_radius', 02172 routineP = moduleN//':'//routineN 02173 02174 CHARACTER(LEN=default_string_length) :: inp_kind_name, kind_name 02175 INTEGER :: handle, i, i_rep, n_rep 02176 LOGICAL :: found 02177 REAL(KIND=dp) :: mm_radius 02178 TYPE(atomic_kind_type), POINTER :: atomic_kind 02179 TYPE(fist_potential_type), POINTER :: fist_potential 02180 TYPE(section_vals_type), POINTER :: kind_section 02181 02182 CALL timeset(routineN,handle) 02183 02184 kind_section => section_vals_get_subs_vals(subsys_section,"KIND",error=error) 02185 CALL section_vals_get(kind_section,n_repetition=n_rep,error=error) 02186 02187 DO i = 1, SIZE(atomic_kind_set) 02188 atomic_kind => atomic_kind_set(i) 02189 CALL get_atomic_kind(atomic_kind=atomic_kind, & 02190 fist_potential=fist_potential, name=kind_name) 02191 CALL uppercase(kind_name) 02192 found = .FALSE. 02193 02194 ! Try to find a matching KIND section in the SUBSYS section and read the 02195 ! MM_RADIUS field if it is present. In case the kind section is never 02196 ! encountered, the mm_radius remains zero. 02197 mm_radius = 0.0_dp 02198 DO i_rep = 1, n_rep 02199 CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_",& 02200 c_val=inp_kind_name, i_rep_section=i_rep, error=error) 02201 CALL uppercase(inp_kind_name) 02202 IF (iw>0) WRITE(iw,*) "Matching kinds for MM_RADIUS :: '", & 02203 TRIM(kind_name), "' with '", TRIM(inp_kind_name), "'" 02204 IF (TRIM(kind_name)==TRIM(inp_kind_name)) THEN 02205 CALL section_vals_val_get(kind_section, i_rep_section=i_rep,& 02206 keyword_name="MM_RADIUS", r_val=mm_radius, error=error) 02207 CALL issue_duplications(found, routineP, "MM_RADIUS", kind_name, error=error) 02208 found = .TRUE. 02209 END IF 02210 END DO 02211 02212 CALL set_potential(potential=fist_potential, mm_radius=mm_radius) 02213 END DO 02214 CALL timestop(handle) 02215 END SUBROUTINE force_field_pack_radius 02216 02217 ! ***************************************************************************** 02221 SUBROUTINE force_field_pack_pol(atomic_kind_set, iw, inp_info, error) 02222 02223 TYPE(atomic_kind_type), DIMENSION(:), 02224 POINTER :: atomic_kind_set 02225 INTEGER, INTENT(IN) :: iw 02226 TYPE(input_info_type), POINTER :: inp_info 02227 TYPE(cp_error_type), INTENT(INOUT) :: error 02228 02229 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_pol', 02230 routineP = moduleN//':'//routineN 02231 02232 CHARACTER(LEN=default_string_length) :: kind_name 02233 INTEGER :: handle, i, j 02234 LOGICAL :: found 02235 REAL(KIND=dp) :: apol, cpol 02236 TYPE(atomic_kind_type), POINTER :: atomic_kind 02237 TYPE(fist_potential_type), POINTER :: fist_potential 02238 02239 CALL timeset(routineN,handle) 02240 02241 DO i=1,SIZE(atomic_kind_set) 02242 atomic_kind => atomic_kind_set(i) 02243 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02244 fist_potential=fist_potential,& 02245 name=kind_name) 02246 CALL get_potential(potential=fist_potential, apol=apol, cpol=cpol) 02247 CALL uppercase(kind_name) 02248 found = .FALSE. 02249 02250 ! Always have the input param last to overwrite all the other ones 02251 IF(ASSOCIATED(inp_info%apol_atm)) THEN 02252 DO j=1,SIZE(inp_info%apol_atm) 02253 IF (iw>0) WRITE(iw,*) "Matching kinds for APOL :: '", & 02254 TRIM(kind_name), "' with '", TRIM(inp_info%apol_atm(j)), "'" 02255 IF((inp_info%apol_atm(j))==kind_name) THEN 02256 apol = inp_info%apol(j) 02257 CALL issue_duplications(found,routineP, "APOL", kind_name, error=error) 02258 found = .TRUE. 02259 END IF 02260 END DO 02261 END IF 02262 02263 IF(ASSOCIATED(inp_info%cpol_atm)) THEN 02264 DO j=1,SIZE(inp_info%cpol_atm) 02265 IF (iw>0) WRITE(iw,*) "Matching kinds for CPOL :: '", & 02266 TRIM(kind_name), "' with '", TRIM(inp_info%cpol_atm(j)), "'" 02267 IF((inp_info%cpol_atm(j))==kind_name) THEN 02268 cpol = inp_info%cpol(j) 02269 CALL issue_duplications(found,routineP, "CPOL", kind_name, error=error) 02270 found = .TRUE. 02271 END IF 02272 END DO 02273 END IF 02274 02275 CALL set_potential(potential=fist_potential, apol=apol, cpol=cpol) 02276 END DO 02277 CALL timestop(handle) 02278 END SUBROUTINE force_field_pack_pol 02279 02280 ! ***************************************************************************** 02283 SUBROUTINE force_field_pack_damp (atomic_kind_set,& 02284 iw,inp_info,error) 02285 02286 TYPE(atomic_kind_type), DIMENSION(:), 02287 POINTER :: atomic_kind_set 02288 INTEGER :: iw 02289 TYPE(input_info_type), POINTER :: inp_info 02290 TYPE(cp_error_type), INTENT(INOUT) :: error 02291 02292 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_damp', 02293 routineP = moduleN//':'//routineN 02294 02295 CHARACTER(len=default_string_length) :: atm_name1, atm_name2, 02296 my_atm_name1, my_atm_name2 02297 INTEGER :: handle2, i, j, k, nkinds 02298 LOGICAL :: found 02299 TYPE(atomic_kind_type), POINTER :: atomic_kind, atomic_kind2 02300 TYPE(damping_p_type), POINTER :: damping 02301 02302 CALL timeset(routineN,handle2) 02303 NULLIFY(damping) 02304 nkinds=SIZE(atomic_kind_set) 02305 02306 #ifdef DEBUG_DAMPING 02307 WRITE(*,*) "DAMPING COUCOU" 02308 ! CALL timestop(handle2) 02309 ! RETURN 02310 #endif 02311 DO j=1,SIZE(atomic_kind_set) 02312 atomic_kind=>atomic_kind_set(j) 02313 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02314 name=atm_name1) 02315 CALL UPPERCASE(atm_name1) 02316 #ifdef DEBUG_DAMPING 02317 WRITE(*,*) "DAMPING in",j,TRIM(atm_name1) 02318 #endif 02319 IF (ASSOCIATED(inp_info%damping_list)) THEN 02320 DO i=1,SIZE(inp_info%damping_list) 02321 my_atm_name1=inp_info%damping_list(i)%atm_name1 02322 my_atm_name2=inp_info%damping_list(i)%atm_name2 02323 02324 IF (iw>0) WRITE(iw,*) "Damping Checking ::",TRIM(my_atm_name1),& 02325 TRIM(atm_name1) 02326 #ifdef DEBUG_DAMPING 02327 WRITE(*,*) "DAMPING DEBUG, checking ::",TRIM(my_atm_name1),& 02328 TRIM(atm_name1) 02329 #endif 02330 IF (my_atm_name1==atm_name1) THEN 02331 IF (.NOT.ASSOCIATED(damping)) THEN 02332 CALL damping_p_create(damping,nkinds,error) 02333 END IF 02334 02335 found=.FALSE. 02336 DO k=1,SIZE(atomic_kind_set) 02337 atomic_kind2=>atomic_kind_set(k) 02338 CALL get_atomic_kind(atomic_kind=atomic_kind2,& 02339 name=atm_name2) 02340 CALL UPPERCASE(atm_name2) 02341 02342 IF (my_atm_name2==atm_name2) THEN 02343 #ifdef DEBUG_DAMPING 02344 WRITE(*,*) "DAMPING FOUND",TRIM(my_atm_name1),& 02345 TRIM(atm_name1),& 02346 TRIM(my_atm_name2),& 02347 TRIM(atm_name2) 02348 #endif 02349 IF(damping%damp(k)%bij/=HUGE(0.0_dp)) found=.TRUE. 02350 CALL issue_duplications(found,routineP, "Damping",& 02351 atm_name1,error=error) 02352 found=.TRUE. 02353 02354 SELECT CASE(TRIM(inp_info%damping_list(i)%dtype)) 02355 CASE('TANG-TOENNIES') 02356 damping%damp(k)%itype=tang_toennies 02357 CASE DEFAULT 02358 CALL stop_program(routineN,moduleN,__LINE__,& 02359 "Unknown damping type.") 02360 END SELECT 02361 damping%damp(k)%order=inp_info%damping_list(i)%order 02362 damping%damp(k)%bij=inp_info%damping_list(i)%bij 02363 damping%damp(k)%cij=inp_info%damping_list(i)%cij 02364 ENDIF 02365 02366 END DO 02367 CALL cp_assert(found,cp_warning_level,cp_assertion_failed,& 02368 routineP,& 02369 "Atom "//TRIM(my_atm_name2)//& 02370 " in damping parameters for atom "//TRIM(my_atm_name1)//& 02371 " not found! "//& 02372 CPSourceFileRef,& 02373 only_ionode=.TRUE.) 02374 ENDIF 02375 END DO 02376 02377 ENDIF 02378 02379 CALL set_atomic_kind(atomic_kind=atomic_kind,damping=damping) 02380 NULLIFY(damping) 02381 END DO 02382 02383 CALL timestop(handle2) 02384 02385 END SUBROUTINE force_field_pack_damp 02386 02387 ! ***************************************************************************** 02390 SUBROUTINE force_field_pack_shell(particle_set, atomic_kind_set,& 02391 molecule_kind_set, molecule_set, root_section, subsys_section,& 02392 shell_particle_set, core_particle_set, cell, iw, inp_info, error) 02393 02394 TYPE(particle_type), DIMENSION(:), 02395 POINTER :: particle_set 02396 TYPE(atomic_kind_type), DIMENSION(:), 02397 POINTER :: atomic_kind_set 02398 TYPE(molecule_kind_type), DIMENSION(:), 02399 POINTER :: molecule_kind_set 02400 TYPE(molecule_type), DIMENSION(:), 02401 POINTER :: molecule_set 02402 TYPE(section_vals_type), POINTER :: root_section, subsys_section 02403 TYPE(particle_type), DIMENSION(:), 02404 POINTER :: shell_particle_set, 02405 core_particle_set 02406 TYPE(cell_type), POINTER :: cell 02407 INTEGER :: iw 02408 TYPE(input_info_type), POINTER :: inp_info 02409 TYPE(cp_error_type), INTENT(INOUT) :: error 02410 02411 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_shell', 02412 routineP = moduleN//':'//routineN 02413 02414 CHARACTER(LEN=default_string_length) :: atmname 02415 INTEGER :: counter, first, first_shell, 02416 handle2, i, j, last, 02417 last_shell, n, natom, nmol, 02418 nshell_tot, stat 02419 INTEGER, DIMENSION(:), POINTER :: molecule_list, shell_list_tmp 02420 LOGICAL :: core_coord_read, failure, found_shell, is_a_shell, 02421 is_link_atom, null_massfrac, only_qm, save_mem, shell_adiabatic, 02422 shell_coord_read 02423 REAL(KIND=dp) :: atmmass 02424 TYPE(atomic_kind_type), POINTER :: atomic_kind 02425 TYPE(molecule_kind_type), POINTER :: molecule_kind 02426 TYPE(molecule_type), POINTER :: molecule 02427 TYPE(section_vals_type), POINTER :: global_section 02428 TYPE(shell_kind_type), POINTER :: shell 02429 TYPE(shell_type), DIMENSION(:), POINTER :: shell_list 02430 02431 failure = .FALSE. 02432 02433 CALL timeset(routineN,handle2) 02434 nshell_tot = 0 02435 n = 0 02436 first_shell = 1 02437 null_massfrac =.FALSE. 02438 core_coord_read = .FALSE. 02439 shell_coord_read = .FALSE. 02440 02441 NULLIFY(global_section) 02442 global_section => section_vals_get_subs_vals(root_section,"GLOBAL",error=error) 02443 CALL section_vals_val_get(global_section,"SAVE_MEM",l_val=save_mem,error=error) 02444 02445 DO i=1,SIZE(atomic_kind_set) 02446 atomic_kind => atomic_kind_set(i) 02447 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02448 name=atmname) 02449 02450 found_shell = .FALSE. 02451 only_qm = qmmm_ff_precond_only_qm(id1=atmname,is_link=is_link_atom) 02452 CALL uppercase(atmname) 02453 02454 ! The shell potential can be defined only from input 02455 IF(ASSOCIATED(inp_info%shell_list)) THEN 02456 DO j=1,SIZE(inp_info%shell_list) 02457 IF (iw > 0) WRITE(iw,*)"Shell Checking ::",TRIM(inp_info%shell_list(j)%atm_name),atmname 02458 IF((inp_info%shell_list(j)%atm_name)==atmname) THEN 02459 CALL get_atomic_kind(atomic_kind=atomic_kind,& 02460 shell=shell, mass=atmmass, natom=natom) 02461 IF(.NOT. ASSOCIATED(shell)) THEN 02462 CALL shell_create(shell,error) 02463 END IF 02464 nshell_tot = nshell_tot + natom 02465 shell%charge_core=inp_info%shell_list(j)%shell%charge_core 02466 shell%charge_shell=inp_info%shell_list(j)%shell%charge_shell 02467 shell%massfrac=inp_info%shell_list(j)%shell%massfrac 02468 IF (shell%massfrac < EPSILON(1.0_dp)) null_massfrac = .TRUE. 02469 shell%k2_spring = inp_info%shell_list(j)%shell%k2_spring 02470 shell%k4_spring = inp_info%shell_list(j)%shell%k4_spring 02471 shell%max_dist=inp_info%shell_list(j)%shell%max_dist 02472 shell%shell_cutoff=inp_info%shell_list(j)%shell%shell_cutoff 02473 shell%mass_shell=shell%massfrac*atmmass 02474 shell%mass_core=atmmass-shell%mass_shell 02475 CALL issue_duplications(found_shell,routineP, "Shell", atmname, error=error) 02476 found_shell = .TRUE. 02477 CALL set_atomic_kind(atomic_kind=atomic_kind,& 02478 shell=shell, shell_active=.TRUE.,error=error) 02479 CALL shell_release(shell, error) 02480 END IF 02481 END DO ! j shell kind 02482 END IF ! associated shell_list 02483 END DO ! i atomic kind 02484 02485 IF (iw > 0) WRITE(iw,*) "Total number of particles with a shell :: ", nshell_tot 02486 ! If shell-model is present: Create particle_set of shells (coord. vel. force) 02487 NULLIFY (shell_particle_set) 02488 NULLIFY (core_particle_set) 02489 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,shell_adiabatic=shell_adiabatic) 02490 IF (nshell_tot > 0) THEN 02491 IF (shell_adiabatic .AND. null_massfrac) THEN 02492 CALL stop_program(routineN,moduleN,__LINE__,& 02493 "Shell-model adiabatic: at least one shell_kind has mass zero") 02494 END IF 02495 CALL allocate_particle_set(shell_particle_set,nshell_tot,error) 02496 CALL allocate_particle_set(core_particle_set,nshell_tot,error) 02497 counter = 0 02498 ! Initialise the shell (and core) coordinates with the particle (atomic) coordinates, 02499 ! count the shell and set pointers 02500 DO i=1,SIZE(particle_set) 02501 NULLIFY (atomic_kind) 02502 NULLIFY (shell) 02503 atomic_kind => particle_set(i)%atomic_kind 02504 CALL get_atomic_kind(atomic_kind=atomic_kind,shell_active=is_a_shell) 02505 IF (is_a_shell) THEN 02506 counter = counter + 1 02507 particle_set(i)%shell_index = counter 02508 shell_particle_set(counter)%shell_index = counter 02509 shell_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind 02510 shell_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3) 02511 shell_particle_set(counter)%atom_index = i 02512 core_particle_set(counter)%shell_index = counter 02513 core_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind 02514 core_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3) 02515 core_particle_set(counter)%atom_index = i 02516 ELSE 02517 particle_set(i)%shell_index = 0 02518 END IF 02519 END DO 02520 CPPostcondition(counter==nshell_tot,cp_failure_level,routineP,error,failure) 02521 END IF 02522 02523 ! Read the shell (and core) coordinates from the restart file, if available 02524 CALL read_binary_cs_coordinates("SHELL",shell_particle_set,root_section,& 02525 subsys_section,shell_coord_read,error) 02526 CALL read_binary_cs_coordinates("CORE",core_particle_set,root_section,& 02527 subsys_section,core_coord_read,error) 02528 02529 IF (nshell_tot > 0) THEN 02530 ! Read the shell (and core) coordinates from the input, if no coordinates were found 02531 ! in the restart file 02532 IF (shell_adiabatic) THEN 02533 IF (.NOT.(core_coord_read.AND.shell_coord_read)) THEN 02534 CALL read_shell_coord_input(particle_set,shell_particle_set,cell,& 02535 subsys_section,core_particle_set,& 02536 save_mem=save_mem,error=error) 02537 END IF 02538 ELSE 02539 IF (.NOT.shell_coord_read) THEN 02540 CALL read_shell_coord_input(particle_set,shell_particle_set,cell,& 02541 subsys_section,save_mem=save_mem,error=error) 02542 END IF 02543 END IF 02544 ! Determine the number of shells per molecule kind 02545 n = 0 02546 DO i=1,SIZE(molecule_kind_set) 02547 molecule_kind => molecule_kind_set(i) 02548 CALL get_molecule_kind(molecule_kind=molecule_kind,molecule_list=molecule_list,& 02549 natom=natom, nmolecule=nmol) 02550 molecule=>molecule_set(molecule_list(1)) 02551 CALL get_molecule(molecule=molecule,first_atom=first,last_atom=last) 02552 ALLOCATE (shell_list_tmp(natom),STAT=stat) 02553 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02554 counter = 0 02555 DO j = first, last 02556 atomic_kind => particle_set(j)%atomic_kind 02557 CALL get_atomic_kind(atomic_kind=atomic_kind,shell_active=is_a_shell) 02558 IF (is_a_shell) THEN 02559 counter = counter+1 02560 shell_list_tmp(counter) = j - first + 1 02561 first_shell = MIN(first_shell, MAX(1,particle_set(j)%shell_index)) 02562 END IF 02563 END DO ! j atom in molecule_kind i, molecule 1 of the molecule_list 02564 IF (counter /= 0) THEN 02565 ! Setup of fist_shell and last_shell for all molecules.. 02566 DO j=1,SIZE(molecule_list) 02567 last_shell = first_shell + counter -1 02568 molecule=>molecule_set(molecule_list(j)) 02569 molecule%first_shell = first_shell 02570 molecule%last_shell = last_shell 02571 first_shell = last_shell + 1 02572 END DO 02573 ! Setup of shell_list 02574 CALL get_molecule_kind(molecule_kind=molecule_kind,shell_list=shell_list) 02575 IF (ASSOCIATED(shell_list)) THEN 02576 DEALLOCATE(shell_list,STAT=stat) 02577 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02578 END IF 02579 ALLOCATE (shell_list(counter), STAT=stat) 02580 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02581 DO j=1,counter 02582 shell_list(j)%a = shell_list_tmp(j) 02583 atomic_kind => particle_set(shell_list_tmp(j)+first-1)%atomic_kind 02584 CALL get_atomic_kind(atomic_kind=atomic_kind,name=atmname, shell=shell) 02585 CALL uppercase(atmname) 02586 shell_list(j)%name = atmname 02587 shell_list(j)%shell_kind => shell 02588 END DO 02589 CALL set_molecule_kind(molecule_kind=molecule_kind, nshell=counter,shell_list=shell_list) 02590 END IF 02591 DEALLOCATE (shell_list_tmp,STAT=stat) 02592 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02593 n = n + nmol*counter 02594 END DO ! i molecule kind 02595 END IF 02596 CPPostcondition(first_shell-1==nshell_tot,cp_failure_level,routineP,error,failure) 02597 CPPostcondition(n==nshell_tot,cp_failure_level,routineP,error,failure) 02598 CALL timestop(handle2) 02599 02600 END SUBROUTINE force_field_pack_shell 02601 02602 ! ***************************************************************************** 02605 SUBROUTINE force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, & 02606 Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env, & 02607 iw2, iw3, error) 02608 02609 TYPE(atomic_kind_type), DIMENSION(:), 02610 POINTER :: atomic_kind_set 02611 TYPE(force_field_type), INTENT(INOUT) :: ff_type 02612 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 02613 INTEGER :: iw 02614 CHARACTER(LEN=default_string_length), 02615 DIMENSION(:), POINTER :: Ainfo 02616 TYPE(charmm_info_type), POINTER :: chm_info 02617 TYPE(input_info_type), POINTER :: inp_info 02618 TYPE(gromos_info_type), POINTER :: gro_info 02619 TYPE(amber_info_type), POINTER :: amb_info 02620 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond14 02621 TYPE(ewald_environment_type), POINTER :: ewald_env 02622 INTEGER :: iw2, iw3 02623 TYPE(cp_error_type), INTENT(INOUT) :: error 02624 02625 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond14', 02626 routineP = moduleN//':'//routineN 02627 02628 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, 02629 name_atm_b, name_atm_b_local 02630 INTEGER :: handle2, i, ii, j, jj, k, 02631 match_names 02632 LOGICAL :: failure, found, found_a, 02633 found_b, only_qm, use_qmmm_ff 02634 REAL(KIND=dp) :: epsilon0, epsilon_a, 02635 epsilon_b, ewald_rcut, rmin, 02636 rmin2_a, rmin2_b 02637 TYPE(atomic_kind_type), POINTER :: atomic_kind 02638 TYPE(pair_potential_single_type), 02639 POINTER :: pot 02640 02641 failure = .FALSE. 02642 use_qmmm_ff = qmmm_env%use_qmmm_ff 02643 NULLIFY(pot) 02644 CALL ewald_env_get(ewald_env, rcut=ewald_rcut, error=error) 02645 CALL timeset(routineN,handle2) 02646 CALL pair_potential_pp_create (potparm_nonbond14, SIZE(atomic_kind_set), error) 02647 DO i=1,SIZE(atomic_kind_set) 02648 atomic_kind => atomic_kind_set(i) 02649 CALL get_atomic_kind(atomic_kind=atomic_kind,name=name_atm_a_local) 02650 DO j=i,SIZE(atomic_kind_set) 02651 atomic_kind => atomic_kind_set(j) 02652 CALL get_atomic_kind(atomic_kind=atomic_kind,name=name_atm_b_local) 02653 found = .FALSE. 02654 found_a = .FALSE. 02655 found_b = .FALSE. 02656 name_atm_a = name_atm_a_local 02657 name_atm_b = name_atm_b_local 02658 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 02659 CALL uppercase(name_atm_a) 02660 CALL uppercase(name_atm_b) 02661 pot => potparm_nonbond14%pot(i,j)%pot 02662 02663 ! loop over params from GROMOS 02664 IF(ASSOCIATED(gro_info%nonbond_a_14)) THEN 02665 ii = 0 02666 jj = 0 02667 DO k=1,SIZE(gro_info%nonbond_a_14) 02668 IF(TRIM(name_atm_a)==TRIM(gro_info%nonbond_a_14(k))) THEN 02669 ii = k 02670 found_a = .TRUE. 02671 EXIT 02672 END IF 02673 END DO 02674 DO k=1,SIZE(gro_info%nonbond_a_14) 02675 IF(TRIM(name_atm_b)==TRIM(gro_info%nonbond_a_14(k))) THEN 02676 jj = k 02677 found_b = .TRUE. 02678 EXIT 02679 END IF 02680 END DO 02681 IF(ii/=0 .AND. jj/=0) THEN 02682 CALL pair_potential_lj_create(pot%set(1)%lj, error) 02683 pot%type = lj_type 02684 pot%at1 = name_atm_a 02685 pot%at2 = name_atm_b 02686 pot%set(1)%lj%epsilon = 1.0_dp 02687 pot%set(1)%lj%sigma6 = gro_info%nonbond_c6_14(ii,jj) 02688 pot%set(1)%lj%sigma12 = gro_info%nonbond_c12_14(ii,jj) 02689 pot%rcutsq = (10.0_dp*bohr)**2 02690 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a, name_atm_b, error=error) 02691 found = .TRUE. 02692 END IF 02693 END IF 02694 02695 ! loop over params from CHARMM 02696 ii = 0 02697 jj = 0 02698 IF(ASSOCIATED(chm_info%nonbond_a_14)) THEN 02699 DO k=1,SIZE(chm_info%nonbond_a_14) 02700 IF((name_atm_a)==(chm_info%nonbond_a_14(k))) THEN 02701 ii = k 02702 rmin2_a = chm_info%nonbond_rmin2_14(k) 02703 epsilon_a = chm_info%nonbond_eps_14(k) 02704 found_a = .TRUE. 02705 END IF 02706 END DO 02707 DO k=1,SIZE(chm_info%nonbond_a_14) 02708 IF((name_atm_b)==(chm_info%nonbond_a_14(k))) THEN 02709 jj = k 02710 rmin2_b = chm_info%nonbond_rmin2_14(k) 02711 epsilon_b = chm_info%nonbond_eps_14(k) 02712 found_b = .TRUE. 02713 END IF 02714 END DO 02715 END IF 02716 IF(ASSOCIATED(chm_info%nonbond_a)) THEN 02717 IF(.NOT.found_a) THEN 02718 DO k=1,SIZE(chm_info%nonbond_a) 02719 IF((name_atm_a)==(chm_info%nonbond_a(k))) THEN 02720 ii = k 02721 rmin2_a = chm_info%nonbond_rmin2(k) 02722 epsilon_a = chm_info%nonbond_eps(k) 02723 END IF 02724 END DO 02725 END IF 02726 IF(.NOT.found_b) THEN 02727 DO k=1,SIZE(chm_info%nonbond_a) 02728 IF((name_atm_b)==(chm_info%nonbond_a(k))) THEN 02729 jj = k 02730 rmin2_b = chm_info%nonbond_rmin2(k) 02731 epsilon_b = chm_info%nonbond_eps(k) 02732 END IF 02733 END DO 02734 END IF 02735 END IF 02736 IF(ii/=0 .AND. jj/=0) THEN 02737 rmin = rmin2_a + rmin2_b 02738 ! ABS to allow for mixing the two different sign conventions for epsilon 02739 epsilon0 = SQRT(ABS(epsilon_a*epsilon_b)) 02740 CALL pair_potential_lj_create(pot%set(1)%lj, error) 02741 pot%type = lj_charmm_type 02742 pot%at1 = name_atm_a 02743 pot%at2 = name_atm_b 02744 pot%set(1)%lj%epsilon = epsilon0 02745 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 02746 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 02747 pot%rcutsq = (10.0_dp*bohr)**2 02748 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a, name_atm_b,& 02749 error=error) 02750 found = .TRUE. 02751 END IF 02752 02753 ! loop over params from AMBER 02754 IF(ASSOCIATED(amb_info%nonbond_a)) THEN 02755 ii = 0 02756 jj = 0 02757 IF(.NOT.found_a) THEN 02758 DO k=1,SIZE(amb_info%nonbond_a) 02759 IF((name_atm_a)==(amb_info%nonbond_a(k))) THEN 02760 ii = k 02761 rmin2_a = amb_info%nonbond_rmin2(k) 02762 epsilon_a = amb_info%nonbond_eps(k) 02763 END IF 02764 END DO 02765 END IF 02766 IF(.NOT.found_b) THEN 02767 DO k=1,SIZE(amb_info%nonbond_a) 02768 IF((name_atm_b)==(amb_info%nonbond_a(k))) THEN 02769 jj = k 02770 rmin2_b = amb_info%nonbond_rmin2(k) 02771 epsilon_b = amb_info%nonbond_eps(k) 02772 END IF 02773 END DO 02774 END IF 02775 IF(ii/=0 .AND. jj/=0) THEN 02776 rmin = rmin2_a + rmin2_b 02777 ! ABS to allow for mixing the two different sign conventions for epsilon 02778 epsilon0 = SQRT(ABS(epsilon_a*epsilon_b)) 02779 CALL pair_potential_lj_create(pot%set(1)%lj, error) 02780 pot%type = lj_charmm_type 02781 pot%at1 = name_atm_a 02782 pot%at2 = name_atm_b 02783 pot%set(1)%lj%epsilon = epsilon0 02784 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 02785 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 02786 pot%rcutsq = (10.0_dp*bohr)**2 02787 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a,& 02788 name_atm_b, error=error) 02789 found = .TRUE. 02790 END IF 02791 END IF 02792 02793 ! always have the input param last to overwrite all the other ones 02794 IF(ASSOCIATED(inp_info%nonbonded14)) THEN 02795 DO k=1,SIZE(inp_info%nonbonded14%pot) 02796 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 02797 " with ",TRIM(inp_info%nonbonded14%pot(k)%pot%at1),& 02798 TRIM(inp_info%nonbonded14%pot(k)%pot%at2) 02799 IF((((name_atm_a)==(inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 02800 ((name_atm_b)==(inp_info%nonbonded14%pot(k)%pot%at2))) .OR. & 02801 (((name_atm_b)==(inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 02802 ((name_atm_a)==(inp_info%nonbonded14%pot(k)%pot%at2))) ) THEN 02803 IF (ff_type%multiple_potential) THEN 02804 CALL pair_potential_single_add(inp_info%nonbonded14%pot(k)%pot,pot,error) 02805 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 02806 "Multiple ONFO declaration: "//TRIM(name_atm_a)//& 02807 " and "//TRIM(name_atm_b)//" ADDING! "//& 02808 CPSourceFileRef,& 02809 only_ionode=.TRUE.) 02810 potparm_nonbond14%pot(i,j)%pot => pot 02811 potparm_nonbond14%pot(j,i)%pot => pot 02812 ELSE 02813 CALL pair_potential_single_copy(inp_info%nonbonded14%pot(k)%pot,pot,error) 02814 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 02815 "Multiple ONFO declarations: "//TRIM(name_atm_a)//& 02816 " and "//TRIM(name_atm_b)//" OVERWRITING! "//& 02817 CPSourceFileRef,& 02818 only_ionode=.TRUE.) 02819 END IF 02820 IF(iw>0) WRITE(iw,*) " FOUND ",TRIM(name_atm_a)," ",TRIM(name_atm_b) 02821 found = .TRUE. 02822 END IF 02823 END DO 02824 END IF 02825 ! at the very end we offer the possibility to overwrite the parameters for QM/MM 02826 ! nonbonded interactions 02827 IF (use_qmmm_ff) THEN 02828 match_names = 0 02829 IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1 02830 IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1 02831 IF (match_names == 1) THEN 02832 IF (ASSOCIATED(qmmm_env%inp_info%nonbonded14)) THEN 02833 DO k=1,SIZE(qmmm_env%inp_info%nonbonded14%pot) 02834 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 02835 " with ",TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1),& 02836 TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2) 02837 IF(( ((name_atm_a) ==(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 02838 ((name_atm_b) ==(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) .OR. & 02839 (((name_atm_b)==(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 02840 ((name_atm_a) ==(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) ) THEN 02841 IF (qmmm_env%multiple_potential) THEN 02842 CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded14%pot(k)%pot,pot,error) 02843 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 02844 "Multiple ONFO declaration: "//TRIM(name_atm_a)//& 02845 " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! "//& 02846 CPSourceFileRef,& 02847 only_ionode=.TRUE.) 02848 potparm_nonbond14%pot(i,j)%pot => pot 02849 potparm_nonbond14%pot(j,i)%pot => pot 02850 ELSE 02851 CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded14%pot(k)%pot,pot,error) 02852 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 02853 "Multiple ONFO declaration: "//TRIM(name_atm_a)//& 02854 " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! "//& 02855 CPSourceFileRef,& 02856 only_ionode=.TRUE.) 02857 END IF 02858 IF(iw>0) WRITE(iw,*) " FOUND ",TRIM(name_atm_a),& 02859 " ",TRIM(name_atm_b) 02860 found = .TRUE. 02861 END IF 02862 END DO 02863 END IF 02864 END IF 02865 END IF 02866 02867 IF(.NOT.found) THEN 02868 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 02869 atm2=TRIM(name_atm_b),& 02870 type_name="Spline_Bond_Env",& 02871 array=Ainfo,& 02872 error=error) 02873 CALL pair_potential_single_clean(pot, error=error) 02874 pot%type = nn_type 02875 pot%at1 = name_atm_a 02876 pot%at2 = name_atm_b 02877 END IF 02878 ! If defined global RCUT let's use it 02879 IF(ff_type%rcut_nb>0.0_dp) THEN 02880 pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb 02881 END IF 02882 ! Cutoff is defined always as the maximum between the FF and Ewald 02883 pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) 02884 IF (only_qm) THEN 02885 CALL pair_potential_single_clean(pot, error=error) 02886 END IF 02887 END DO ! atom kind j 02888 END DO ! atom kind i 02889 CALL timestop(handle2) 02890 02891 END SUBROUTINE force_field_pack_nonbond14 02892 02893 ! ***************************************************************************** 02896 SUBROUTINE force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, fatal, & 02897 iw, Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond, & 02898 ewald_env, error) 02899 02900 TYPE(atomic_kind_type), DIMENSION(:), 02901 POINTER :: atomic_kind_set 02902 TYPE(force_field_type), INTENT(INOUT) :: ff_type 02903 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 02904 LOGICAL :: fatal 02905 INTEGER :: iw 02906 CHARACTER(LEN=default_string_length), 02907 DIMENSION(:), POINTER :: Ainfo 02908 TYPE(charmm_info_type), POINTER :: chm_info 02909 TYPE(input_info_type), POINTER :: inp_info 02910 TYPE(gromos_info_type), POINTER :: gro_info 02911 TYPE(amber_info_type), POINTER :: amb_info 02912 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond 02913 TYPE(ewald_environment_type), POINTER :: ewald_env 02914 TYPE(cp_error_type), INTENT(INOUT) :: error 02915 02916 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond', 02917 routineP = moduleN//':'//routineN 02918 02919 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, 02920 name_atm_b, name_atm_b_local 02921 INTEGER :: handle2, i, ii, j, jj, k, 02922 match_names 02923 LOGICAL :: failure, found, is_a_shell, 02924 is_b_shell, only_qm, 02925 use_qmmm_ff 02926 REAL(KIND=dp) :: epsilon0, ewald_rcut, rmin 02927 TYPE(atomic_kind_type), POINTER :: atomic_kind 02928 TYPE(pair_potential_single_type), 02929 POINTER :: pot 02930 02931 CALL timeset(routineN,handle2) 02932 failure = .FALSE. 02933 use_qmmm_ff = qmmm_env%use_qmmm_ff 02934 NULLIFY(pot) 02935 CALL ewald_env_get(ewald_env, rcut=ewald_rcut, error=error) 02936 CALL pair_potential_pp_create ( potparm_nonbond, SIZE(atomic_kind_set),error ) 02937 DO i=1,SIZE(atomic_kind_set) 02938 atomic_kind => atomic_kind_set(i) 02939 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local, & 02940 shell_active=is_a_shell) 02941 DO j=i,SIZE(atomic_kind_set) 02942 atomic_kind => atomic_kind_set(j) 02943 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local, & 02944 shell_active=is_b_shell) 02945 found = .FALSE. 02946 name_atm_a = name_atm_a_local 02947 name_atm_b = name_atm_b_local 02948 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 02949 CALL uppercase(name_atm_a) 02950 CALL uppercase(name_atm_b) 02951 pot => potparm_nonbond%pot(i,j)%pot 02952 02953 ! loop over params from GROMOS 02954 IF(ASSOCIATED(gro_info%nonbond_a)) THEN 02955 ii = 0 02956 jj = 0 02957 DO k=1,SIZE(gro_info%nonbond_a) 02958 IF(TRIM(name_atm_a)==TRIM(gro_info%nonbond_a(k))) THEN 02959 ii = k 02960 EXIT 02961 END IF 02962 END DO 02963 DO k=1,SIZE(gro_info%nonbond_a) 02964 IF(TRIM(name_atm_b)==TRIM(gro_info%nonbond_a(k))) THEN 02965 jj = k 02966 EXIT 02967 END IF 02968 END DO 02969 02970 IF(ii/=0 .AND. jj/=0) THEN 02971 CALL pair_potential_lj_create(pot%set(1)%lj, error) 02972 pot%type = lj_type 02973 pot%at1 = name_atm_a 02974 pot%at2 = name_atm_b 02975 pot%set(1)%lj%epsilon = 1.0_dp 02976 pot%set(1)%lj%sigma6 = gro_info%nonbond_c6(ii,jj) 02977 pot%set(1)%lj%sigma12 = gro_info%nonbond_c12(ii,jj) 02978 pot%rcutsq = (10.0_dp*bohr)**2 02979 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a, name_atm_b, error=error) 02980 found = .TRUE. 02981 END IF 02982 END IF 02983 02984 ! loop over params from CHARMM 02985 IF(ASSOCIATED(chm_info%nonbond_a)) THEN 02986 ii = 0 02987 jj = 0 02988 DO k=1,SIZE(chm_info%nonbond_a) 02989 IF((name_atm_a)==(chm_info%nonbond_a(k))) THEN 02990 ii = k 02991 END IF 02992 END DO 02993 DO k=1,SIZE(chm_info%nonbond_a) 02994 IF((name_atm_b)==(chm_info%nonbond_a(k))) THEN 02995 jj = k 02996 END IF 02997 END DO 02998 02999 IF(ii/=0 .AND. jj/=0) THEN 03000 rmin = chm_info%nonbond_rmin2(ii)+chm_info%nonbond_rmin2(jj) 03001 epsilon0 = SQRT(chm_info%nonbond_eps(ii)*& 03002 chm_info%nonbond_eps(jj)) 03003 CALL pair_potential_lj_create(pot%set(1)%lj, error) 03004 pot%type = lj_charmm_type 03005 pot%at1 = name_atm_a 03006 pot%at2 = name_atm_b 03007 pot%set(1)%lj%epsilon = epsilon0 03008 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 03009 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 03010 pot%rcutsq = (10.0_dp*bohr)**2 03011 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a, name_atm_b, error=error) 03012 found = .TRUE. 03013 END IF 03014 END IF 03015 03016 ! loop over params from AMBER 03017 IF(ASSOCIATED(amb_info%nonbond_a)) THEN 03018 ii = 0 03019 jj = 0 03020 DO k=1,SIZE(amb_info%nonbond_a) 03021 IF((name_atm_a)==(amb_info%nonbond_a(k))) THEN 03022 ii = k 03023 END IF 03024 END DO 03025 DO k=1,SIZE(amb_info%nonbond_a) 03026 IF((name_atm_b)==(amb_info%nonbond_a(k))) THEN 03027 jj = k 03028 END IF 03029 END DO 03030 03031 IF(ii/=0 .AND. jj/=0) THEN 03032 rmin = amb_info%nonbond_rmin2(ii)+amb_info%nonbond_rmin2(jj) 03033 epsilon0 = SQRT(amb_info%nonbond_eps(ii)*amb_info%nonbond_eps(jj)) 03034 CALL pair_potential_lj_create(pot%set(1)%lj, error) 03035 pot%type = lj_charmm_type 03036 pot%at1 = name_atm_a 03037 pot%at2 = name_atm_b 03038 pot%set(1)%lj%epsilon = epsilon0 03039 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 03040 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 03041 pot%rcutsq = (10.0_dp*bohr)**2 03042 CALL issue_duplications(found,routineP, "Lennard-Jones", name_atm_a, name_atm_b, error=error) 03043 found = .TRUE. 03044 END IF 03045 END IF 03046 03047 ! always have the input param last to overwrite all the other ones 03048 IF(ASSOCIATED(inp_info%nonbonded)) THEN 03049 DO k=1,SIZE(inp_info%nonbonded%pot) 03050 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1)=="*").OR.& 03051 (TRIM(inp_info%nonbonded%pot(k)%pot%at2)=="*")) CYCLE 03052 03053 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 03054 " with ",TRIM(inp_info%nonbonded%pot(k)%pot%at1),& 03055 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 03056 IF(( ((name_atm_a)==(inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 03057 ((name_atm_b)==(inp_info%nonbonded%pot(k)%pot%at2))) .OR. & 03058 (((name_atm_b)==(inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 03059 ((name_atm_a)==(inp_info%nonbonded%pot(k)%pot%at2))) ) THEN 03060 IF (ff_type%multiple_potential) THEN 03061 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot,pot,error) 03062 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03063 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03064 " and "//TRIM(name_atm_b)//" ADDING! "//& 03065 CPSourceFileRef,& 03066 only_ionode=.TRUE.) 03067 potparm_nonbond%pot(i,j)%pot => pot 03068 potparm_nonbond%pot(j,i)%pot => pot 03069 ELSE 03070 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot,pot,error) 03071 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03072 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03073 " and "//TRIM(name_atm_b)//" OVERWRITING! "//& 03074 CPSourceFileRef,& 03075 only_ionode=.TRUE.) 03076 END IF 03077 IF(iw>0) WRITE(iw,*) " FOUND ",TRIM(name_atm_a)," ",TRIM(name_atm_b) 03078 found = .TRUE. 03079 END IF 03080 END DO 03081 ! Check for wildcards for one of the two types (if not associated yet) 03082 IF (.not.found) THEN 03083 DO k=1,SIZE(inp_info%nonbonded%pot) 03084 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1)=="*").EQV.& 03085 (TRIM(inp_info%nonbonded%pot(k)%pot%at2)=="*")) CYCLE 03086 03087 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 03088 " with ",TRIM(inp_info%nonbonded%pot(k)%pot%at1),& 03089 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 03090 03091 IF( (name_atm_a==inp_info%nonbonded%pot(k)%pot%at1) .OR.& 03092 (name_atm_b==inp_info%nonbonded%pot(k)%pot%at2) .OR.& 03093 (name_atm_b==inp_info%nonbonded%pot(k)%pot%at1) .OR.& 03094 (name_atm_a==inp_info%nonbonded%pot(k)%pot%at2) ) THEN 03095 IF (ff_type%multiple_potential) THEN 03096 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot,pot,error) 03097 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03098 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03099 " and "//TRIM(name_atm_b)//" ADDING! "//& 03100 CPSourceFileRef,& 03101 only_ionode=.TRUE.) 03102 potparm_nonbond%pot(i,j)%pot => pot 03103 potparm_nonbond%pot(j,i)%pot => pot 03104 ELSE 03105 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot,pot,error) 03106 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03107 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03108 " and "//TRIM(name_atm_b)//" OVERWRITING! "//& 03109 CPSourceFileRef,& 03110 only_ionode=.TRUE.) 03111 END IF 03112 IF(iw>0) WRITE(iw,*) " FOUND (one WILDCARD)",TRIM(name_atm_a)," ",TRIM(name_atm_b) 03113 found = .TRUE. 03114 END IF 03115 END DO 03116 END IF 03117 ! Check for wildcards for both types (if not associated yet) 03118 IF (.not.found) THEN 03119 DO k=1,SIZE(inp_info%nonbonded%pot) 03120 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1)/="*").OR.& 03121 (TRIM(inp_info%nonbonded%pot(k)%pot%at2)/="*")) CYCLE 03122 03123 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 03124 " with ",TRIM(inp_info%nonbonded%pot(k)%pot%at1),& 03125 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 03126 03127 IF (ff_type%multiple_potential) THEN 03128 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot,pot,error) 03129 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03130 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03131 " and "//TRIM(name_atm_b)//" ADDING! "//& 03132 CPSourceFileRef,& 03133 only_ionode=.TRUE.) 03134 potparm_nonbond%pot(i,j)%pot => pot 03135 potparm_nonbond%pot(j,i)%pot => pot 03136 ELSE 03137 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot,pot,error) 03138 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03139 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03140 " and "//TRIM(name_atm_b)//" OVERWRITING! "//& 03141 CPSourceFileRef,& 03142 only_ionode=.TRUE.) 03143 END IF 03144 IF(iw>0) WRITE(iw,*) " FOUND (both WILDCARDS)",TRIM(name_atm_a)," ",TRIM(name_atm_b) 03145 found = .TRUE. 03146 END DO 03147 END IF 03148 END IF 03149 03150 ! at the very end we offer the possibility to overwrite the parameters for QM/MM 03151 ! nonbonded interactions 03152 IF (use_qmmm_ff) THEN 03153 match_names = 0 03154 IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1 03155 IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1 03156 IF (match_names == 1) THEN 03157 IF (ASSOCIATED(qmmm_env%inp_info%nonbonded)) THEN 03158 DO k=1,SIZE(qmmm_env%inp_info%nonbonded%pot) 03159 IF(iw>0) WRITE(iw,*) " TESTING ",TRIM(name_atm_a),TRIM(name_atm_b),& 03160 " with ",TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1),& 03161 TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2) 03162 IF(( ((name_atm_a) ==(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 03163 ((name_atm_b) ==(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) .OR. & 03164 (((name_atm_b)==(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 03165 ((name_atm_a) ==(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) ) THEN 03166 IF (qmmm_env%multiple_potential) THEN 03167 CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded%pot(k)%pot,pot,error) 03168 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03169 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03170 " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! "//& 03171 CPSourceFileRef,& 03172 only_ionode=.TRUE.) 03173 potparm_nonbond%pot(i,j)%pot => pot 03174 potparm_nonbond%pot(j,i)%pot => pot 03175 ELSE 03176 CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded%pot(k)%pot,pot,error) 03177 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routineP,& 03178 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)//& 03179 " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! "//& 03180 CPSourceFileRef,& 03181 only_ionode=.TRUE.) 03182 END IF 03183 IF(iw>0) WRITE(iw,*) " FOUND ",TRIM(name_atm_a)," ",TRIM(name_atm_b) 03184 found = .TRUE. 03185 END IF 03186 END DO 03187 END IF 03188 END IF 03189 END IF 03190 IF(.NOT.found) THEN 03191 CALL store_FF_missing_par(atm1=TRIM(name_atm_a),& 03192 atm2=TRIM(name_atm_b),& 03193 type_name="Spline_Non_Bond_Env",& 03194 fatal=fatal,& 03195 array=Ainfo,& 03196 error=error) 03197 END IF 03198 ! If defined global RCUT let's use it 03199 IF(ff_type%rcut_nb>0.0_dp) THEN 03200 pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb 03201 END IF 03202 ! Cutoff is defined always as the maximum between the FF and Ewald 03203 pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) 03204 ! Set the shell type 03205 IF((is_a_shell .AND. .NOT. is_b_shell) .OR. (is_b_shell .AND. .NOT. is_a_shell)) THEN 03206 pot%shell_type = nosh_sh 03207 ELSE IF (is_a_shell .AND. is_b_shell) THEN 03208 pot%shell_type = sh_sh 03209 ELSE 03210 pot%shell_type = nosh_nosh 03211 END IF 03212 IF (only_qm) THEN 03213 CALL pair_potential_single_clean(pot, error=error) 03214 END IF 03215 END DO ! jkind 03216 END DO ! ikind 03217 CALL timestop(handle2) 03218 END SUBROUTINE force_field_pack_nonbond 03219 03220 ! ***************************************************************************** 03223 SUBROUTINE force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, & 03224 potparm, do_zbl, nonbonded_type, error) 03225 03226 TYPE(atomic_kind_type), DIMENSION(:), 03227 POINTER :: atomic_kind_set 03228 TYPE(force_field_type), INTENT(INOUT) :: ff_type 03229 INTEGER :: iw2, iw3, iw4 03230 TYPE(pair_potential_pp_type), POINTER :: potparm 03231 LOGICAL, INTENT(IN) :: do_zbl 03232 CHARACTER(LEN=*), INTENT(IN) :: nonbonded_type 03233 TYPE(cp_error_type), INTENT(INOUT) :: error 03234 03235 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_splines', 03236 routineP = moduleN//':'//routineN 03237 03238 INTEGER :: handle2, ikind, jkind, n 03239 LOGICAL :: failure 03240 TYPE(spline_data_p_type), DIMENSION(:), 03241 POINTER :: spl_p 03242 TYPE(spline_environment_type), POINTER :: spline_env 03243 03244 CALL timeset(routineN,handle2) 03245 failure = .FALSE. 03246 ! Figure out which nonbonded interactions happen to be indentical, and 03247 ! prepare storage for these, avoiding duplicates. 03248 NULLIFY(spline_env) 03249 CALL get_nonbond_storage(spline_env, potparm, atomic_kind_set, & 03250 do_zbl, shift_cutoff=ff_type%shift_cutoff, error=error) 03251 ! Effectively compute the spline data. 03252 CALL spline_nonbond_control(spline_env, potparm, & 03253 atomic_kind_set, eps_spline=ff_type%eps_spline, & 03254 max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb, & 03255 emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, iw=iw2, iw2=iw3, iw3=iw4, & 03256 do_zbl=do_zbl, shift_cutoff=ff_type%shift_cutoff, & 03257 nonbonded_type=nonbonded_type, error=error) 03258 ! Let the pointers on potparm point to the splines generated in 03259 ! spline_nonbond_control. 03260 DO ikind = 1, SIZE ( potparm%pot, 1 ) 03261 DO jkind = ikind, SIZE ( potparm%pot, 2) 03262 n = spline_env % spltab ( ikind, jkind ) 03263 spl_p => spline_env%spl_pp(n)%spl_p 03264 CALL spline_data_p_retain ( spl_p, error ) 03265 CALL spline_data_p_release ( potparm%pot(ikind,jkind)%pot%pair_spline_data, error ) 03266 potparm%pot(ikind,jkind)%pot%pair_spline_data => spl_p 03267 END DO 03268 END DO 03269 CALL spline_env_release(spline_env,error) 03270 CALL timestop(handle2) 03271 03272 END SUBROUTINE force_field_pack_splines 03273 03274 ! ***************************************************************************** 03278 SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, & 03279 potparm_nonbond, ewald_env, error) 03280 03281 TYPE(atomic_kind_type), DIMENSION(:), 03282 POINTER :: atomic_kind_set 03283 TYPE(force_field_type), INTENT(IN) :: ff_type 03284 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond 03285 TYPE(ewald_environment_type), POINTER :: ewald_env 03286 TYPE(cp_error_type), INTENT(INOUT) :: error 03287 03288 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_eicut', 03289 routineP = moduleN//':'//routineN 03290 03291 INTEGER :: ewald_type, handle, i1, i2, 03292 nkinds, stat 03293 LOGICAL :: failure 03294 REAL(KIND=dp) :: alpha, beta, mm_radius1, 03295 mm_radius2, rcut2, 03296 rcut2_ewald, tmp 03297 REAL(KIND=dp), DIMENSION(:, :, :), 03298 POINTER :: interaction_cutoffs 03299 TYPE(atomic_kind_type), POINTER :: atomic_kind 03300 03301 CALL timeset(routineN,handle) 03302 failure = .FALSE. 03303 03304 tmp = 0.0_dp 03305 nkinds = SIZE(atomic_kind_set) 03306 ! allocate the array with interaction cutoffs for the electrostatics, used 03307 ! to make the electrostatic interaction continuous at ewald_env%rcut 03308 ALLOCATE(interaction_cutoffs(3,nkinds,nkinds), stat=stat) 03309 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03310 interaction_cutoffs = 0.0_dp 03311 03312 ! compute the interaction cutoff if SHIFT_CUTOFF is active 03313 IF (ff_type%shift_cutoff) THEN 03314 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, & 03315 rcut=rcut2_ewald) 03316 rcut2_ewald = rcut2_ewald*rcut2_ewald 03317 DO i1 = 1, nkinds 03318 atomic_kind => atomic_kind_set(i1) 03319 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius1) 03320 DO i2 = 1, nkinds 03321 rcut2=rcut2_ewald 03322 IF (ASSOCIATED(potparm_nonbond)) THEN 03323 rcut2 = MAX(potparm_nonbond%pot(i1,i2)%pot%rcutsq, rcut2_ewald) 03324 END IF 03325 IF (rcut2 > 0) THEN 03326 atomic_kind => atomic_kind_set(i2) 03327 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius2) 03328 ! cutoff for core-core 03329 interaction_cutoffs(1,i1,i2) = potential_coulomb(rcut2, tmp, & 03330 1.0_dp, ewald_type, alpha, 0.0_dp, 0.0_dp) 03331 ! cutoff for core-shell, core-ion, shell-core or ion-core 03332 IF (mm_radius1 > 0.0_dp) THEN 03333 beta = sqrthalf/mm_radius1 03334 ELSE 03335 beta = 0.0_dp 03336 END IF 03337 interaction_cutoffs(2,i1,i2) = potential_coulomb(rcut2, tmp, & 03338 1.0_dp, ewald_type, alpha, beta, 0.0_dp) 03339 ! cutoff for shell-shell or ion-ion 03340 IF (mm_radius1 + mm_radius2 > 0.0_dp) THEN 03341 beta = sqrthalf/SQRT(mm_radius1*mm_radius1+mm_radius2*mm_radius2) 03342 ELSE 03343 beta = 0.0_dp 03344 END IF 03345 interaction_cutoffs(3,i1,i2) = potential_coulomb(rcut2, tmp, & 03346 1.0_dp, ewald_type, alpha, beta, 0.0_dp) 03347 END IF 03348 END DO 03349 END DO 03350 END IF 03351 03352 CALL ewald_env_set(ewald_env, interaction_cutoffs=interaction_cutoffs) 03353 03354 CALL timestop(handle) 03355 END SUBROUTINE force_field_pack_eicut 03356 03357 ! ***************************************************************************** 03362 SUBROUTINE issue_duplications(found, routinePext, tag_label, name_atm_a, name_atm_b,& 03363 name_atm_c, name_atm_d, error) 03364 03365 LOGICAL, INTENT(IN) :: found 03366 CHARACTER(LEN=*), INTENT(IN) :: routinePext, tag_label, 03367 name_atm_a 03368 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name_atm_b, name_atm_c, 03369 name_atm_d 03370 TYPE(cp_error_type), INTENT(INOUT) :: error 03371 03372 CHARACTER(len=*), PARAMETER :: routineN = 'issue_duplications', 03373 routineP = moduleN//':'//routineN 03374 03375 CHARACTER(LEN=default_string_length) :: item 03376 LOGICAL :: failure 03377 03378 failure = .FALSE. 03379 item="( "//TRIM(name_atm_a) 03380 IF (PRESENT(name_atm_b)) THEN 03381 item = TRIM(item)//" , "//TRIM(name_atm_b) 03382 END IF 03383 IF (PRESENT(name_atm_c)) THEN 03384 item = TRIM(item)//" , "//TRIM(name_atm_c) 03385 END IF 03386 IF (PRESENT(name_atm_d)) THEN 03387 item = TRIM(item)//" , "//TRIM(name_atm_d) 03388 END IF 03389 item = TRIM(item)//" )" 03390 CALL cp_assert(.NOT.found,cp_warning_level,cp_assertion_failed,routinePext,& 03391 "Multiple "//TRIM(tag_label)//" declarations: "//TRIM(item)//" overwriting! "//& 03392 CPSourceFileRef,& 03393 only_ionode=.TRUE.) 03394 03395 END SUBROUTINE issue_duplications 03396 03397 ! ***************************************************************************** 03400 SUBROUTINE store_FF_missing_par(atm1,atm2,atm3,atm4,type_name,fatal,array,error) 03401 CHARACTER(LEN=*), INTENT(IN) :: atm1 03402 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: atm2, atm3, atm4 03403 CHARACTER(LEN=*), INTENT(IN) :: type_name 03404 LOGICAL, INTENT(INOUT), OPTIONAL :: fatal 03405 CHARACTER(LEN=default_string_length), 03406 DIMENSION(:), POINTER :: array 03407 TYPE(cp_error_type), INTENT(inout) :: error 03408 03409 CHARACTER(len=*), PARAMETER :: routineN = 'store_FF_missing_par', 03410 routineP = moduleN//':'//routineN 03411 03412 CHARACTER(LEN=10) :: sfmt 03413 CHARACTER(LEN=4) :: my_atm1, my_atm2, my_atm3, 03414 my_atm4 03415 CHARACTER(LEN=default_path_length) :: my_format 03416 INTEGER :: fmt, i, nsize 03417 LOGICAL :: failure, found 03418 03419 failure = .FALSE. 03420 nsize = 0 03421 fmt = 1 03422 my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)//& 03423 '",T40,"(",A4,")")' 03424 IF (PRESENT(atm2)) fmt = fmt+1 03425 IF (PRESENT(atm3)) fmt = fmt+1 03426 IF (PRESENT(atm4)) fmt = fmt+1 03427 CALL integer_to_string(fmt-1,sfmt) 03428 IF (fmt>1) & 03429 my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)//& 03430 '",T40,"(",A4,'//TRIM(sfmt)//'(",",A4),")")' 03431 IF (.NOT.failure) THEN 03432 IF (PRESENT(fatal)) fatal = .TRUE. 03433 ! Check for previous already stored equal force fields 03434 IF (ASSOCIATED(array)) nsize = SIZE(array) 03435 found = .FALSE. 03436 IF (nsize>=1) THEN 03437 DO i = 1, nsize 03438 SELECT CASE (type_name) 03439 CASE("Bond") 03440 IF (INDEX(array(i)(21:39),"Bond")==0) CYCLE 03441 my_atm1 = array(i)(41:44) 03442 my_atm2 = array(i)(46:49) 03443 CALL compress(my_atm1,.TRUE.) 03444 CALL compress(my_atm2,.TRUE.) 03445 IF (((atm1==my_atm1).AND.(atm2==my_atm2)).OR.& 03446 ((atm1==my_atm2).AND.(atm2==my_atm1))) found = .TRUE. 03447 CASE("Angle") 03448 IF (INDEX(array(i)(21:39),"Angle")==0) CYCLE 03449 my_atm1 = array(i)(41:44) 03450 my_atm2 = array(i)(46:49) 03451 my_atm3 = array(i)(51:54) 03452 CALL compress(my_atm1,.TRUE.) 03453 CALL compress(my_atm2,.TRUE.) 03454 CALL compress(my_atm3,.TRUE.) 03455 IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3)).OR.& 03456 ((atm1==my_atm3).AND.(atm2==my_atm2).AND.(atm3==my_atm1)))& 03457 found = .TRUE. 03458 CASE("Urey-Bradley") 03459 IF (INDEX(array(i)(21:39),"Urey-Bradley")==0) CYCLE 03460 my_atm1 = array(i)(41:44) 03461 my_atm2 = array(i)(46:49) 03462 my_atm3 = array(i)(51:54) 03463 CALL compress(my_atm1,.TRUE.) 03464 CALL compress(my_atm2,.TRUE.) 03465 CALL compress(my_atm3,.TRUE.) 03466 IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3)).OR.& 03467 ((atm1==my_atm3).AND.(atm2==my_atm2).AND.(atm3==my_atm1)))& 03468 found = .TRUE. 03469 CASE("Torsion") 03470 IF (INDEX(array(i)(21:39),"Torsion")==0) CYCLE 03471 my_atm1 = array(i)(41:44) 03472 my_atm2 = array(i)(46:49) 03473 my_atm3 = array(i)(51:54) 03474 my_atm4 = array(i)(56:59) 03475 CALL compress(my_atm1,.TRUE.) 03476 CALL compress(my_atm2,.TRUE.) 03477 CALL compress(my_atm3,.TRUE.) 03478 CALL compress(my_atm4,.TRUE.) 03479 IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3).AND.(atm4==my_atm4)).OR.& 03480 ((atm1==my_atm4).AND.(atm2==my_atm3).AND.(atm3==my_atm2).AND.(atm4==my_atm1)))& 03481 found = .TRUE. 03482 CASE("Improper") 03483 IF (INDEX(array(i)(21:39),"Improper")==0) CYCLE 03484 my_atm1 = array(i)(41:44) 03485 my_atm2 = array(i)(46:49) 03486 my_atm3 = array(i)(51:54) 03487 my_atm4 = array(i)(56:59) 03488 CALL compress(my_atm1,.TRUE.) 03489 CALL compress(my_atm2,.TRUE.) 03490 CALL compress(my_atm3,.TRUE.) 03491 CALL compress(my_atm4,.TRUE.) 03492 IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3).AND.(atm4==my_atm4)).OR.& 03493 ((atm1==my_atm1).AND.(atm2==my_atm3).AND.(atm3==my_atm2).AND.(atm4==my_atm4)).OR.& 03494 ((atm1==my_atm1).AND.(atm2==my_atm3).AND.(atm3==my_atm4).AND.(atm4==my_atm3)).OR.& 03495 ((atm1==my_atm1).AND.(atm2==my_atm4).AND.(atm3==my_atm3).AND.(atm4==my_atm2)).OR.& 03496 ((atm1==my_atm1).AND.(atm2==my_atm4).AND.(atm3==my_atm2).AND.(atm4==my_atm3)).OR.& 03497 ((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm4).AND.(atm4==my_atm3)))& 03498 found = .TRUE. 03499 03500 CASE("Out of plane bend") 03501 IF (INDEX(array(i)(21:39),"Out of plane bend")==0) CYCLE 03502 my_atm1 = array(i)(41:44) 03503 my_atm2 = array(i)(46:49) 03504 my_atm3 = array(i)(51:54) 03505 my_atm4 = array(i)(56:59) 03506 CALL compress(my_atm1,.TRUE.) 03507 CALL compress(my_atm2,.TRUE.) 03508 CALL compress(my_atm3,.TRUE.) 03509 CALL compress(my_atm4,.TRUE.) 03510 IF (((atm1==my_atm1).AND.(atm2==my_atm2).AND.(atm3==my_atm3).AND.(atm4==my_atm4)).OR.& 03511 ((atm1==my_atm1).AND.(atm2==my_atm3).AND.(atm3==my_atm2).AND.(atm4==my_atm4)))& 03512 found = .TRUE. 03513 03514 CASE("Charge") 03515 IF (INDEX(array(i)(21:39),"Charge")==0) CYCLE 03516 my_atm1 = array(i)(41:44) 03517 CALL compress(my_atm1,.TRUE.) 03518 IF (atm1==my_atm1) found = .TRUE. 03519 CASE("Spline_Bond_Env","Spline_Non_Bond_Env") 03520 IF (INDEX(array(i)(21:39),"Spline_")==0) CYCLE 03521 fmt = 0 03522 my_atm1 = array(i)(41:44) 03523 my_atm2 = array(i)(46:49) 03524 CALL compress(my_atm1,.TRUE.) 03525 CALL compress(my_atm2,.TRUE.) 03526 IF (((atm1==my_atm1).AND.(atm2==my_atm2)).OR.& 03527 ((atm1==my_atm2).AND.(atm2==my_atm1))) found = .TRUE. 03528 CASE DEFAULT 03529 ! Should never reach this point 03530 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 03531 END SELECT 03532 IF (found) EXIT 03533 END DO 03534 ENDIF 03535 IF (.NOT.found) THEN 03536 nsize = nsize + 1 03537 CALL reallocate(array,1,nsize) 03538 SELECT CASE(fmt) 03539 CASE(1) 03540 WRITE(array(nsize),FMT=TRIM(my_format))atm1 03541 CASE(2) 03542 WRITE(array(nsize),FMT=TRIM(my_format))atm1,atm2 03543 CASE(3) 03544 WRITE(array(nsize),FMT=TRIM(my_format))atm1,atm2,atm3 03545 CASE(4) 03546 WRITE(array(nsize),FMT=TRIM(my_format))atm1,atm2,atm3,atm4 03547 END SELECT 03548 END IF 03549 END IF 03550 03551 END SUBROUTINE store_FF_missing_par 03552 03553 END MODULE force_fields_all 03554
1.7.3