1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Define the data structure for the molecule information. 8!> \par History 9!> JGH (22.05.2004) add last_atom information 10!> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints 11!> (patch by Marcel Baer) 12!> \author MK (29.08.2003) 13! ************************************************************************************************** 14MODULE molecule_types 15 16 USE colvar_types, ONLY: colvar_counters,& 17 colvar_release,& 18 colvar_type 19 USE kinds, ONLY: dp 20 USE molecule_kind_types, ONLY: colvar_constraint_type,& 21 fixd_constraint_type,& 22 g3x3_constraint_type,& 23 g4x6_constraint_type,& 24 get_molecule_kind,& 25 molecule_kind_type,& 26 vsite_constraint_type 27#include "../base/base_uses.f90" 28 29 IMPLICIT NONE 30 31 PRIVATE 32 33! *** Global parameters (in this module) *** 34 35 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_types' 36 37! *** Data types *** 38! ************************************************************************************************** 39 TYPE local_colvar_constraint_type 40 LOGICAL :: init 41 TYPE(colvar_type), POINTER :: colvar 42 TYPE(colvar_type), POINTER :: colvar_old 43 REAL(KIND=dp) :: lambda, sigma 44 END TYPE local_colvar_constraint_type 45 46! ************************************************************************************************** 47 TYPE local_g3x3_constraint_type 48 LOGICAL :: init 49 REAL(KIND=dp) :: scale, scale_old, imass1, imass2, imass3 50 REAL(KIND=dp), DIMENSION(3) :: fa, fb, fc, f_roll1, f_roll2, f_roll3, & 51 ra_old, rb_old, rc_old, & 52 va, vb, vc, lambda, del_lambda, lambda_old, & 53 r0_12, r0_13, r0_23 54 REAL(KIND=dp), DIMENSION(3, 3) :: amat 55 END TYPE local_g3x3_constraint_type 56 57! ************************************************************************************************** 58 TYPE local_g4x6_constraint_type 59 LOGICAL :: init 60 REAL(KIND=dp) :: scale, scale_old, imass1, imass2, imass3, imass4 61 REAL(KIND=dp), DIMENSION(3) :: fa, fb, fc, fd, fe, ff, & 62 f_roll1, f_roll2, f_roll3, f_roll4, f_roll5, f_roll6, & 63 ra_old, rb_old, rc_old, rd_old, re_old, rf_old, & 64 va, vb, vc, vd, ve, vf, & 65 r0_12, r0_13, r0_14, r0_23, r0_24, r0_34 66 REAL(KIND=dp), DIMENSION(6) :: lambda, del_lambda, lambda_old 67 REAL(KIND=dp), DIMENSION(6, 6) :: amat 68 END TYPE local_g4x6_constraint_type 69 70! ************************************************************************************************** 71 TYPE local_states_type 72 INTEGER, DIMENSION(:), POINTER :: states ! indices of Kohn-Sham states for molecule 73 INTEGER :: nstates ! Kohn-Sham states for molecule 74 END TYPE local_states_type 75 76! ************************************************************************************************** 77 TYPE local_constraint_type 78 TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv 79 TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 80 TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 81 END TYPE local_constraint_type 82 83! ************************************************************************************************** 84 TYPE global_constraint_type 85 TYPE(colvar_counters) :: ncolv 86 INTEGER :: ntot, nrestraint 87 INTEGER :: ng3x3, ng3x3_restraint 88 INTEGER :: ng4x6, ng4x6_restraint 89 INTEGER :: nvsite, nvsite_restraint 90 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 91 TYPE(colvar_constraint_type), DIMENSION(:), POINTER :: colv_list 92 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list 93 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list 94 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list 95 TYPE(local_colvar_constraint_type), DIMENSION(:), POINTER :: lcolv 96 TYPE(local_g3x3_constraint_type), DIMENSION(:), POINTER :: lg3x3 97 TYPE(local_g4x6_constraint_type), DIMENSION(:), POINTER :: lg4x6 98 END TYPE global_constraint_type 99 100! ************************************************************************************************** 101 TYPE molecule_type 102 TYPE(molecule_kind_type), POINTER :: molecule_kind ! pointer to molecule kind information 103 TYPE(local_states_type), DIMENSION(:), POINTER :: lmi ! local (spin)-states information 104 TYPE(local_constraint_type), POINTER :: lci ! local molecule constraint info 105 INTEGER :: first_atom ! global index of first atom in molecule 106 INTEGER :: last_atom ! global index of last atom in molecule 107 INTEGER :: first_shell ! global index of first shell atom in molecule 108 INTEGER :: last_shell ! global index of last shell atom in molecule 109 END TYPE molecule_type 110 111! *** Public data types *** 112 113 PUBLIC :: local_colvar_constraint_type, & 114 local_g3x3_constraint_type, & 115 local_g4x6_constraint_type, & 116 local_constraint_type, & 117 local_states_type, & 118 global_constraint_type, & 119 molecule_type 120 121! *** Public subroutines *** 122 123 PUBLIC :: deallocate_global_constraint, & 124 allocate_molecule_set, & 125 deallocate_molecule_set, & 126 get_molecule, & 127 set_molecule, & 128 set_molecule_set, & 129 molecule_of_atom, & 130 get_molecule_set_info 131 132CONTAINS 133 134! ************************************************************************************************** 135!> \brief Deallocate a global constraint. 136!> \param gci ... 137!> \par History 138!> 07.2003 created [fawzi] 139!> 01.2014 moved from cp_subsys_release() into separate routine. 140!> \author Ole Schuett 141! ************************************************************************************************** 142 SUBROUTINE deallocate_global_constraint(gci) 143 TYPE(global_constraint_type), POINTER :: gci 144 145 INTEGER :: i 146 147 IF (ASSOCIATED(gci)) THEN 148 ! List of constraints 149 IF (ASSOCIATED(gci%colv_list)) THEN 150 DO i = 1, SIZE(gci%colv_list) 151 DEALLOCATE (gci%colv_list(i)%i_atoms) 152 END DO 153 DEALLOCATE (gci%colv_list) 154 END IF 155 156 IF (ASSOCIATED(gci%g3x3_list)) & 157 DEALLOCATE (gci%g3x3_list) 158 159 IF (ASSOCIATED(gci%g4x6_list)) & 160 DEALLOCATE (gci%g4x6_list) 161 162 ! Local information 163 IF (ASSOCIATED(gci%lcolv)) THEN 164 DO i = 1, SIZE(gci%lcolv) 165 CALL colvar_release(gci%lcolv(i)%colvar) 166 CALL colvar_release(gci%lcolv(i)%colvar_old) 167 NULLIFY (gci%lcolv(i)%colvar, gci%lcolv(i)%colvar_old) 168 END DO 169 DEALLOCATE (gci%lcolv) 170 ENDIF 171 172 IF (ASSOCIATED(gci%lg3x3)) & 173 DEALLOCATE (gci%lg3x3) 174 175 IF (ASSOCIATED(gci%lg4x6)) & 176 DEALLOCATE (gci%lg4x6) 177 178 IF (ASSOCIATED(gci%fixd_list)) & 179 DEALLOCATE (gci%fixd_list) 180 181 DEALLOCATE (gci) 182 ENDIF 183 END SUBROUTINE deallocate_global_constraint 184 185! ************************************************************************************************** 186!> \brief Allocate a molecule set. 187!> \param molecule_set ... 188!> \param nmolecule ... 189!> \date 29.08.2003 190!> \author MK 191!> \version 1.0 192! ************************************************************************************************** 193 SUBROUTINE allocate_molecule_set(molecule_set, nmolecule) 194 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 195 INTEGER, INTENT(IN) :: nmolecule 196 197 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_molecule_set', & 198 routineP = moduleN//':'//routineN 199 200 INTEGER :: imolecule 201 202 IF (ASSOCIATED(molecule_set)) CALL deallocate_molecule_set(molecule_set) 203 204 ALLOCATE (molecule_set(nmolecule)) 205 206 DO imolecule = 1, nmolecule 207 NULLIFY (molecule_set(imolecule)%molecule_kind) 208 NULLIFY (molecule_set(imolecule)%lmi) 209 NULLIFY (molecule_set(imolecule)%lci) 210 211 molecule_set(imolecule)%first_atom = 0 212 molecule_set(imolecule)%last_atom = 0 213 molecule_set(imolecule)%first_shell = 0 214 molecule_set(imolecule)%last_shell = 0 215 END DO 216 217 END SUBROUTINE allocate_molecule_set 218 219! ************************************************************************************************** 220!> \brief Deallocate a molecule set. 221!> \param molecule_set ... 222!> \date 29.08.2003 223!> \author MK 224!> \version 1.0 225! ************************************************************************************************** 226 SUBROUTINE deallocate_molecule_set(molecule_set) 227 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 228 229 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_molecule_set', & 230 routineP = moduleN//':'//routineN 231 232 INTEGER :: imolecule, j 233 234 IF (ASSOCIATED(molecule_set)) THEN 235 236 DO imolecule = 1, SIZE(molecule_set) 237 IF (ASSOCIATED(molecule_set(imolecule)%lmi)) THEN 238 DO j = 1, SIZE(molecule_set(imolecule)%lmi) 239 IF (ASSOCIATED(molecule_set(imolecule)%lmi(j)%states)) THEN 240 DEALLOCATE (molecule_set(imolecule)%lmi(j)%states) 241 ENDIF 242 END DO 243 DEALLOCATE (molecule_set(imolecule)%lmi) 244 ENDIF 245 IF (ASSOCIATED(molecule_set(imolecule)%lci)) THEN 246 IF (ASSOCIATED(molecule_set(imolecule)%lci%lcolv)) THEN 247 DO j = 1, SIZE(molecule_set(imolecule)%lci%lcolv) 248 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar) 249 CALL colvar_release(molecule_set(imolecule)%lci%lcolv(j)%colvar_old) 250 NULLIFY (molecule_set(imolecule)%lci%lcolv(j)%colvar) 251 NULLIFY (molecule_set(imolecule)%lci%lcolv(j)%colvar_old) 252 END DO 253 DEALLOCATE (molecule_set(imolecule)%lci%lcolv) 254 ENDIF 255 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg3x3)) THEN 256 DEALLOCATE (molecule_set(imolecule)%lci%lg3x3) 257 ENDIF 258 IF (ASSOCIATED(molecule_set(imolecule)%lci%lg4x6)) THEN 259 DEALLOCATE (molecule_set(imolecule)%lci%lg4x6) 260 ENDIF 261 DEALLOCATE (molecule_set(imolecule)%lci) 262 ENDIF 263 ENDDO 264 DEALLOCATE (molecule_set) 265 266 ELSE 267 268 CALL cp_abort(__LOCATION__, & 269 "The pointer molecule_set is not associated and "// & 270 "cannot be deallocated") 271 272 END IF 273 274 END SUBROUTINE deallocate_molecule_set 275 276! ************************************************************************************************** 277!> \brief Get components from a molecule data set. 278!> \param molecule ... 279!> \param molecule_kind ... 280!> \param lmi ... 281!> \param lci ... 282!> \param lg3x3 ... 283!> \param lg4x6 ... 284!> \param lcolv ... 285!> \param first_atom ... 286!> \param last_atom ... 287!> \param first_shell ... 288!> \param last_shell ... 289!> \date 29.08.2003 290!> \author MK 291!> \version 1.0 292! ************************************************************************************************** 293 SUBROUTINE get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, & 294 first_atom, last_atom, first_shell, last_shell) 295 296 TYPE(molecule_type), POINTER :: molecule 297 TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind 298 TYPE(local_states_type), DIMENSION(:), OPTIONAL, & 299 POINTER :: lmi 300 TYPE(local_constraint_type), OPTIONAL, POINTER :: lci 301 TYPE(local_g3x3_constraint_type), OPTIONAL, & 302 POINTER :: lg3x3(:) 303 TYPE(local_g4x6_constraint_type), OPTIONAL, & 304 POINTER :: lg4x6(:) 305 TYPE(local_colvar_constraint_type), DIMENSION(:), & 306 OPTIONAL, POINTER :: lcolv 307 INTEGER, OPTIONAL :: first_atom, last_atom, first_shell, & 308 last_shell 309 310 CHARACTER(len=*), PARAMETER :: routineN = 'get_molecule', routineP = moduleN//':'//routineN 311 312 IF (ASSOCIATED(molecule)) THEN 313 314 IF (PRESENT(first_atom)) first_atom = molecule%first_atom 315 IF (PRESENT(last_atom)) last_atom = molecule%last_atom 316 IF (PRESENT(first_shell)) first_shell = molecule%first_shell 317 IF (PRESENT(last_shell)) last_shell = molecule%last_shell 318 IF (PRESENT(molecule_kind)) molecule_kind => molecule%molecule_kind 319 IF (PRESENT(lmi)) lmi => molecule%lmi 320 IF (PRESENT(lci)) lci => molecule%lci 321 IF (PRESENT(lcolv)) THEN 322 IF (ASSOCIATED(molecule%lci)) THEN 323 lcolv => molecule%lci%lcolv 324 ELSE 325 CPABORT("The pointer lci is not associated") 326 ENDIF 327 ENDIF 328 IF (PRESENT(lg3x3)) THEN 329 IF (ASSOCIATED(molecule%lci)) THEN 330 lg3x3 => molecule%lci%lg3x3 331 ELSE 332 CPABORT("The pointer lci is not associated") 333 ENDIF 334 ENDIF 335 IF (PRESENT(lg4x6)) THEN 336 IF (ASSOCIATED(molecule%lci)) THEN 337 lg4x6 => molecule%lci%lg4x6 338 ELSE 339 CPABORT("The pointer lci is not associated") 340 ENDIF 341 ENDIF 342 343 ELSE 344 345 CPABORT("The pointer lci is not associated") 346 347 END IF 348 349 END SUBROUTINE get_molecule 350 351! ************************************************************************************************** 352!> \brief Set a molecule data set. 353!> \param molecule ... 354!> \param molecule_kind ... 355!> \param lmi ... 356!> \param lci ... 357!> \param lcolv ... 358!> \param lg3x3 ... 359!> \param lg4x6 ... 360!> \date 29.08.2003 361!> \author MK 362!> \version 1.0 363! ************************************************************************************************** 364 SUBROUTINE set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6) 365 TYPE(molecule_type), POINTER :: molecule 366 TYPE(molecule_kind_type), OPTIONAL, POINTER :: molecule_kind 367 TYPE(local_states_type), DIMENSION(:), OPTIONAL, & 368 POINTER :: lmi 369 TYPE(local_constraint_type), OPTIONAL, POINTER :: lci 370 TYPE(local_colvar_constraint_type), DIMENSION(:), & 371 OPTIONAL, POINTER :: lcolv 372 TYPE(local_g3x3_constraint_type), OPTIONAL, & 373 POINTER :: lg3x3(:) 374 TYPE(local_g4x6_constraint_type), OPTIONAL, & 375 POINTER :: lg4x6(:) 376 377 CHARACTER(len=*), PARAMETER :: routineN = 'set_molecule', routineP = moduleN//':'//routineN 378 379 IF (ASSOCIATED(molecule)) THEN 380 381 IF (PRESENT(molecule_kind)) molecule%molecule_kind => molecule_kind 382 IF (PRESENT(lmi)) molecule%lmi => lmi 383 IF (PRESENT(lci)) molecule%lci => lci 384 IF (PRESENT(lcolv)) THEN 385 IF (ASSOCIATED(molecule%lci)) THEN 386 molecule%lci%lcolv => lcolv 387 ELSE 388 CPABORT("The pointer lci is not associated") 389 ENDIF 390 ENDIF 391 IF (PRESENT(lg3x3)) THEN 392 IF (ASSOCIATED(molecule%lci)) THEN 393 molecule%lci%lg3x3 => lg3x3 394 ELSE 395 CPABORT("The pointer lci is not associated") 396 ENDIF 397 ENDIF 398 IF (PRESENT(lg4x6)) THEN 399 IF (ASSOCIATED(molecule%lci)) THEN 400 molecule%lci%lg4x6 => lg4x6 401 ELSE 402 CPABORT("The pointer lci is not associated") 403 ENDIF 404 ENDIF 405 ELSE 406 407 CPABORT("The pointer molecule is not associated") 408 409 END IF 410 411 END SUBROUTINE set_molecule 412 413! ************************************************************************************************** 414!> \brief Set a molecule data set. 415!> \param molecule_set ... 416!> \param first_atom ... 417!> \param last_atom ... 418!> \date 29.08.2003 419!> \author MK 420!> \version 1.0 421! ************************************************************************************************** 422 SUBROUTINE set_molecule_set(molecule_set, first_atom, last_atom) 423 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 424 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: first_atom, last_atom 425 426 CHARACTER(len=*), PARAMETER :: routineN = 'set_molecule_set', & 427 routineP = moduleN//':'//routineN 428 429 INTEGER :: imolecule 430 431 IF (ASSOCIATED(molecule_set)) THEN 432 433 IF (PRESENT(first_atom)) THEN 434 435 IF (SIZE(first_atom) /= SIZE(molecule_set)) THEN 436 CALL cp_abort(__LOCATION__, & 437 "The sizes of first_atom and molecule_set "// & 438 "are different") 439 END IF 440 441 DO imolecule = 1, SIZE(molecule_set) 442 molecule_set(imolecule)%first_atom = first_atom(imolecule) 443 END DO 444 445 END IF 446 447 IF (PRESENT(last_atom)) THEN 448 449 IF (SIZE(last_atom) /= SIZE(molecule_set)) THEN 450 CALL cp_abort(__LOCATION__, & 451 "The sizes of last_atom and molecule_set "// & 452 "are different") 453 END IF 454 455 DO imolecule = 1, SIZE(molecule_set) 456 molecule_set(imolecule)%last_atom = last_atom(imolecule) 457 END DO 458 459 END IF 460 461 ELSE 462 463 CPABORT("The pointer molecule_set is not associated") 464 465 END IF 466 467 END SUBROUTINE set_molecule_set 468 469! ************************************************************************************************** 470!> \brief finds for each atom the molecule it belongs to 471!> \param molecule_set ... 472!> \param atom_to_mol ... 473! ************************************************************************************************** 474 SUBROUTINE molecule_of_atom(molecule_set, atom_to_mol) 475 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 476 INTEGER, DIMENSION(:), INTENT(OUT) :: atom_to_mol 477 478 CHARACTER(len=*), PARAMETER :: routineN = 'molecule_of_atom', & 479 routineP = moduleN//':'//routineN 480 481 INTEGER :: first_atom, iatom, imol, last_atom 482 TYPE(molecule_type), POINTER :: molecule 483 484 DO imol = 1, SIZE(molecule_set) 485 molecule => molecule_set(imol) 486 CALL get_molecule(molecule=molecule, first_atom=first_atom, last_atom=last_atom) 487 DO iatom = first_atom, last_atom 488 atom_to_mol(iatom) = imol 489 ENDDO ! iatom 490 END DO ! imol 491 492 END SUBROUTINE molecule_of_atom 493 494! ************************************************************************************************** 495!> \brief returns information about molecules in the set. 496!> \param molecule_set ... 497!> \param atom_to_mol ... 498!> \param mol_to_first_atom ... 499!> \param mol_to_last_atom ... 500!> \param mol_to_nelectrons ... 501!> \param mol_to_nbasis ... 502!> \param mol_to_charge ... 503!> \param mol_to_multiplicity ... 504!> \par History 505!> 2011.06 created [Rustam Z Khaliullin] 506!> \author Rustam Z Khaliullin 507! ************************************************************************************************** 508 SUBROUTINE get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, & 509 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, & 510 mol_to_multiplicity) 511 512 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 513 INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: atom_to_mol, mol_to_first_atom, & 514 mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity 515 516 CHARACTER(len=*), PARAMETER :: routineN = 'get_molecule_set_info', & 517 routineP = moduleN//':'//routineN 518 519 INTEGER :: first_atom, iatom, imol, last_atom, & 520 nbasis, nelec 521 REAL(KIND=dp) :: charge 522 TYPE(molecule_kind_type), POINTER :: imol_kind 523 TYPE(molecule_type), POINTER :: molecule 524 525 DO imol = 1, SIZE(molecule_set) 526 527 molecule => molecule_set(imol) 528 CALL get_molecule(molecule=molecule, molecule_kind=imol_kind, & 529 first_atom=first_atom, last_atom=last_atom) 530 531 IF (PRESENT(mol_to_nelectrons)) THEN 532 CALL get_molecule_kind(imol_kind, nelectron=nelec) 533 mol_to_nelectrons(imol) = nelec 534 ENDIF 535 536 IF (PRESENT(mol_to_multiplicity)) THEN 537 ! RZK-warning: At the moment we can only get the total number 538 ! of electrons (alpha+beta) and we do not have a way to get the multiplicity of mols. 539 ! Therefore, the best we can do is to assume the singlet state for even number of electrons 540 ! and doublet state for odd number of electrons (assume ne_alpha > ne_beta). 541 ! The best way to implement a correct multiplicity subroutine in the future is to get 542 ! the number of alpha and beta e- for each atom from init_atom_electronic_state. This way (as opposed to 543 ! reading the multiplicities from file) the number of occupied and virtual orbitals 544 ! will be consistent with atomic guess. A guess with broken symmetry will be easy to 545 ! implement as well. 546 CALL get_molecule_kind(imol_kind, nelectron=nelec) 547 IF (MOD(nelec, 2) .EQ. 0) THEN 548 mol_to_multiplicity(imol) = 1 549 ELSE 550 mol_to_multiplicity(imol) = 2 551 ENDIF 552 ENDIF 553 554 IF (PRESENT(mol_to_charge)) THEN 555 CALL get_molecule_kind(imol_kind, charge=charge) 556 mol_to_charge(imol) = NINT(charge) 557 ENDIF 558 559 IF (PRESENT(mol_to_nbasis)) THEN 560 CALL get_molecule_kind(imol_kind, nsgf=nbasis) 561 mol_to_nbasis(imol) = nbasis 562 ENDIF 563 564 IF (PRESENT(mol_to_first_atom)) THEN 565 mol_to_first_atom(imol) = first_atom 566 ENDIF 567 568 IF (PRESENT(mol_to_last_atom)) THEN 569 mol_to_last_atom(imol) = last_atom 570 ENDIF 571 572 IF (PRESENT(atom_to_mol)) THEN 573 DO iatom = first_atom, last_atom 574 atom_to_mol(iatom) = imol 575 ENDDO ! iatom 576 ENDIF 577 578 END DO ! imol 579 580 END SUBROUTINE get_molecule_set_info 581 582END MODULE molecule_types 583