61 USE util,
ONLY: find_boundary,&
63 #include "./base/base_uses.f90"
67 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'topology_generate_util'
70 LOGICAL,
PARAMETER :: debug_this_module = .false.
96 INTEGER,
INTENT(IN) :: natom, natom_prev, nbond_prev
97 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: id_molname
99 CHARACTER(LEN=default_string_length),
PARAMETER :: basename =
"MOL"
101 CHARACTER(LEN=default_string_length) :: molname
102 INTEGER :: i, n, nmol
109 ALLOCATE (atom_bond_list(natom))
111 ALLOCATE (atom_bond_list(i)%array1(0))
114 IF (
ASSOCIATED(conn_info%bond_a)) n =
SIZE(conn_info%bond_a) - nbond_prev
115 CALL reorder_structure(atom_bond_list, conn_info%bond_a(nbond_prev + 1:) - natom_prev, &
116 conn_info%bond_b(nbond_prev + 1:) - natom_prev, n)
119 check = all(id_molname ==
str2id(
s2s(
"__UNDEF__"))) .OR. all(id_molname /=
str2id(
s2s(
"__UNDEF__")))
122 IF (
id2str(id_molname(i)) ==
"__UNDEF__")
THEN
123 molname = trim(basename)//adjustl(cp_to_string(nmol))
124 CALL generate_molname_low(i, atom_bond_list, molname, id_molname)
129 DEALLOCATE (atom_bond_list(i)%array1)
131 DEALLOCATE (atom_bond_list)
144 RECURSIVE SUBROUTINE generate_molname_low(i, atom_bond_list, molname, id_molname)
145 INTEGER,
INTENT(IN) :: i
147 CHARACTER(LEN=default_string_length),
INTENT(IN) :: molname
148 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: id_molname
152 IF (debug_this_module)
THEN
153 WRITE (*, *)
"Entered with :", i
154 WRITE (*, *) trim(molname)//
": entering with i:", i,
" full series to test:: ", atom_bond_list(i)%array1
155 IF ((trim(
id2str(id_molname(i))) /=
"__UNDEF__") .AND. &
156 (trim(
id2str(id_molname(i))) /= trim(molname)))
THEN
157 WRITE (*, *)
"Atom (", i,
") has already a molecular name assigned ! ("//trim(
id2str(id_molname(i)))//
")."
158 WRITE (*, *)
"New molecular name would be: ("//trim(molname)//
")."
159 WRITE (*, *)
"Detecting something wrong in the molecular setup!"
163 id_molname(i) =
str2id(molname)
164 DO j = 1,
SIZE(atom_bond_list(i)%array1)
165 k = atom_bond_list(i)%array1(j)
166 IF (debug_this_module)
WRITE (*, *)
"entering with i:", i,
"testing :", k
168 atom_bond_list(i)%array1(j) = -1
169 WHERE (atom_bond_list(k)%array1 == i) atom_bond_list(k)%array1 = -1
170 CALL generate_molname_low(k, atom_bond_list, molname, id_molname)
172 END SUBROUTINE generate_molname_low
183 LOGICAL,
INTENT(in),
OPTIONAL :: qmmm
184 TYPE(qmmm_env_mm_type),
OPTIONAL,
POINTER :: qmmm_env
185 TYPE(section_vals_type),
POINTER :: subsys_section
187 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_molecule'
188 INTEGER,
PARAMETER :: nblock = 100
190 INTEGER :: atom_in_kind, atom_in_mol, first, handle, handle2, i, iatm, iatom, iend, ifirst, &
191 ilast, inum, istart, itype, iw, j, jump1, jump2, last, max_mol_num, mol_num, mol_res, &
192 mol_typ, myind, n, natom, nlocl, ntype, resid
193 INTEGER,
DIMENSION(:),
POINTER :: qm_atom_index, wrk1, wrk2
194 LOGICAL :: do_again, found, my_qmmm
198 TYPE(cp_logger_type),
POINTER :: logger
203 extension=
".subsysLog")
204 CALL timeset(routinen, handle)
205 NULLIFY (qm_atom_index)
215 IF (
PRESENT(qmmm) .AND.
PRESENT(qmmm_env)) my_qmmm = qmmm
218 IF (
ASSOCIATED(atom_info%map_mol_typ))
DEALLOCATE (atom_info%map_mol_typ)
219 ALLOCATE (atom_info%map_mol_typ(natom))
221 IF (
ASSOCIATED(atom_info%map_mol_num))
DEALLOCATE (atom_info%map_mol_num)
222 ALLOCATE (atom_info%map_mol_num(natom))
224 IF (
ASSOCIATED(atom_info%map_mol_res))
DEALLOCATE (atom_info%map_mol_res)
225 ALLOCATE (atom_info%map_mol_res(natom))
228 atom_info%map_mol_typ(:) = 0
229 atom_info%map_mol_num(:) = -1
230 atom_info%map_mol_res(:) = 1
234 atom_info%map_mol_typ(1) = 1
236 CALL reallocate(wrk1, 1, nblock)
237 wrk1(1) = atom_info%id_molname(1)
242 atom_info%map_mol_typ(iatom) = ntype
246 IF (atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1))
THEN
247 atom_info%map_mol_typ(iatom) = atom_info%map_mol_typ(iatom - 1)
248 IF (atom_info%id_resname(iatom) == atom_info%id_resname(iatom - 1))
THEN
249 atom_info%map_mol_res(iatom) = atom_info%map_mol_res(iatom - 1)
252 atom_info%map_mol_res(iatom) = resid
258 IF (atom_info%id_molname(iatom) == wrk1(itype))
THEN
259 atom_info%map_mol_typ(iatom) = itype
264 IF (.NOT. found)
THEN
266 atom_info%map_mol_typ(iatom) = ntype
267 IF (ntype >
SIZE(wrk1))
CALL reallocate(wrk1, 1, 2*
SIZE(wrk1))
268 wrk1(ntype) = atom_info%id_molname(iatom)
271 atom_info%map_mol_res(iatom) = resid
274 IF (atom_info%id_molname(iatom - 1) == atom_info%id_molname(iatom))
THEN
275 atom_info%map_mol_typ(iatom) = ntype
278 atom_info%map_mol_typ(iatom) = ntype
284 IF (iw > 0)
WRITE (iw,
'(/,T2,A)')
"Start of molecule generation"
288 ALLOCATE (atom_bond_list(natom))
290 ALLOCATE (atom_bond_list(i)%array1(0))
293 IF (
ASSOCIATED(conn_info%bond_a)) n =
SIZE(conn_info%bond_a)
294 CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, n)
295 CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
297 DEALLOCATE (atom_bond_list(i)%array1)
299 DEALLOCATE (atom_bond_list)
300 IF (iw > 0)
WRITE (iw,
'(/,T2,A)')
"End of molecule generation"
303 IF (iw > 0)
WRITE (iw,
'(/,T2,A)')
"Checking for non-continuous generated molecules"
305 ALLOCATE (wrk1(natom))
306 ALLOCATE (wrk2(natom))
307 wrk1 = atom_info%map_mol_num
309 IF (debug_this_module)
THEN
311 WRITE (*,
'(2I10)') i, atom_info%map_mol_num(i)
315 CALL sort(wrk1, natom, wrk2)
317 mol_typ = wrk1(istart)
319 IF (mol_typ /= wrk1(i))
THEN
321 first = minval(wrk2(istart:iend))
322 last = maxval(wrk2(istart:iend))
323 nlocl = last - first + 1
324 IF (iend - istart + 1 /= nlocl)
THEN
325 IF (debug_this_module)
WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl
326 CALL cp_abort(__location__, &
327 "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// &
328 "In particular a molecule defined from index ("//cp_to_string(first)//
") to ("// &
329 cp_to_string(last)//
") contains other molecules, not connected! "// &
330 "Too late at this stage everything should be already ordered! "// &
331 "If you have not yet employed the REORDERING keyword, please do so. "// &
332 "It may help to fix this issue.")
335 mol_typ = wrk1(istart)
339 first = minval(wrk2(istart:iend))
340 last = maxval(wrk2(istart:iend))
341 nlocl = last - first + 1
342 IF (iend - istart + 1 /= nlocl)
THEN
343 IF (debug_this_module)
WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl
344 CALL cp_abort(__location__, &
345 "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// &
346 "In particular a molecule defined from index ("//cp_to_string(first)//
") to ("// &
347 cp_to_string(last)//
") contains other molecules, not connected! "// &
348 "Too late at this stage everything should be already ordered! "// &
349 "If you have not yet employed the REORDERING keyword, please do so. "// &
350 "It may help to fix this issue.")
354 IF (iw > 0)
WRITE (iw,
'(/,T2,A)')
"End of check"
356 IF (iw > 0)
WRITE (unit=iw, fmt=
"(/,T2,A)")
"Start of renumbering molecules"
359 atom_info%map_mol_num(1) = 1
361 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1))
THEN
363 ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1))
THEN
364 mol_num = mol_num + 1
366 atom_info%map_mol_num(iatom) = mol_num
369 mol_typ = atom_info%map_mol_typ(1)
370 mol_num = atom_info%map_mol_num(1)
372 IF (atom_info%map_mol_typ(i) /= mol_typ)
THEN
373 myind = atom_info%map_mol_num(i) - mol_num + 1
374 cpassert(myind /= atom_info%map_mol_num(i - 1))
375 mol_typ = atom_info%map_mol_typ(i)
376 mol_num = atom_info%map_mol_num(i)
378 atom_info%map_mol_num(i) = atom_info%map_mol_num(i) - mol_num + 1
381 IF (iw > 0)
WRITE (unit=iw, fmt=
"(/,T2,A)")
"End of renumbering molecules"
384 CALL timeset(routinen//
"_PARA_RES", handle2)
385 IF (iw > 0)
WRITE (unit=iw, fmt=
"(/,T2,A,L2)")
"Starting PARA_RES: ",
topology%para_res
388 atom_info%id_molname(:) = atom_info%id_resname(:)
390 atom_info%map_mol_typ(1) = 1
392 atom_info%map_mol_num(1) = 1
394 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1))
THEN
397 ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1))
THEN
398 mol_num = mol_num + 1
400 atom_info%map_mol_typ(iatom) = ntype
401 atom_info%map_mol_num(iatom) = mol_num
405 mol_typ = atom_info%map_mol_typ(1)
406 mol_num = atom_info%map_mol_num(1)
407 atom_info%map_mol_res(1) = mol_res
409 IF ((atom_info%resid(i - 1) /= atom_info%resid(i)) .OR. &
410 (atom_info%id_resname(i - 1) /= atom_info%id_resname(i)))
THEN
411 mol_res = mol_res + 1
413 IF ((atom_info%map_mol_typ(i) /= mol_typ) .OR. &
414 (atom_info%map_mol_num(i) /= mol_num))
THEN
415 mol_typ = atom_info%map_mol_typ(i)
416 mol_num = atom_info%map_mol_num(i)
419 atom_info%map_mol_res(i) = mol_res
423 IF (iw > 0)
WRITE (unit=iw, fmt=
"(/,T2,A)")
"End of PARA_RES"
424 CALL timestop(handle2)
428 WRITE (iw,
'(4(1X,A,":",I0),2(1X,A,1X,A))')
"iatom", iatom, &
429 "map_mol_typ", atom_info%map_mol_typ(iatom), &
430 "map_mol_num", atom_info%map_mol_num(iatom), &
431 "map_mol_res", atom_info%map_mol_res(iatom), &
432 "mol_name:", trim(
id2str(atom_info%id_molname(iatom))), &
433 "res_name:", trim(
id2str(atom_info%id_resname(iatom)))
439 IF (iw > 0)
WRITE (iw, *)
"MAP_MOL_NUM ", atom_info%map_mol_num
440 IF (iw > 0)
WRITE (iw, *)
"MAP_MOL_TYP ", atom_info%map_mol_typ
441 IF (iw > 0)
WRITE (iw, *)
"MAP_MOL_RES ", atom_info%map_mol_res
442 ALLOCATE (qm_atom_index(
SIZE(qmmm_env%qm_atom_index)))
443 qm_atom_index = qmmm_env%qm_atom_index
444 cpassert(all(qm_atom_index /= 0))
445 DO myind = 1,
SIZE(qm_atom_index)
446 IF (qm_atom_index(myind) == 0) cycle
447 CALL find_boundary(atom_info%map_mol_typ, natom, ifirst, ilast, &
448 atom_info%map_mol_typ(qm_atom_index(myind)))
449 CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, &
450 atom_info%map_mol_typ(qm_atom_index(myind)), atom_info%map_mol_num(qm_atom_index(myind)))
451 IF (iw > 0)
WRITE (iw, *)
"qm fragment:: ifirst, ilast", ifirst, ilast
452 cpassert(((ifirst /= 0) .OR. (ilast /= natom)))
453 DO iatm = ifirst, ilast
454 atom_info%id_molname(iatm) =
str2id(
s2s(
"_QM_"// &
455 trim(
id2str(atom_info%id_molname(iatm)))))
456 IF (iw > 0)
WRITE (iw, *)
"QM Molecule name :: ",
id2str(atom_info%id_molname(iatm))
457 WHERE (qm_atom_index == iatm) qm_atom_index = 0
459 DO iatm = 1, ifirst - 1
460 IF (any(qm_atom_index == iatm)) do_again = .true.
462 DO iatm = ilast + 1, natom
463 IF (any(qm_atom_index == iatm)) do_again = .true.
465 IF (iw > 0)
WRITE (iw, *)
" Another QM fragment? :: ", do_again
466 IF (ifirst /= 1)
THEN
467 jump1 = atom_info%map_mol_typ(ifirst) - atom_info%map_mol_typ(ifirst - 1)
468 cpassert(jump1 <= 1 .AND. jump1 >= 0)
469 jump1 = abs(jump1 - 1)
473 IF (ilast /= natom)
THEN
474 jump2 = atom_info%map_mol_typ(ilast + 1) - atom_info%map_mol_typ(ilast)
475 cpassert(jump2 <= 1 .AND. jump2 >= 0)
476 jump2 = abs(jump2 - 1)
482 DO iatm = ifirst, natom
483 atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump1
485 DO iatm = ilast + 1, natom
486 atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump2
489 DO iatm = ifirst, ilast
490 atom_info%map_mol_num(iatm) = 1
495 CALL find_boundary(atom_info%map_mol_typ, natom, first, last, atom_info%map_mol_typ(ilast + 1))
496 CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, &
497 atom_info%map_mol_typ(ilast + 1), atom_info%map_mol_num(ilast + 1))
498 atom_in_mol = ilast - ifirst + 1
500 DO iatm = first, last, atom_in_mol
501 atom_info%map_mol_num(iatm:iatm + atom_in_mol - 1) = inum
506 IF (.NOT. do_again)
EXIT
508 DEALLOCATE (qm_atom_index)
511 WRITE (iw, *)
"After the QM/MM Setup:"
513 WRITE (iw, *)
" iatom,map_mol_typ,map_mol_num ", iatom, &
514 atom_info%map_mol_typ(iatom), atom_info%map_mol_num(iatom)
522 WRITE (iw, *)
"SUMMARY:: Number of molecule kinds found:", ntype
523 ntype = maxval(atom_info%map_mol_typ)
525 atom_in_kind = count(atom_info%map_mol_typ == i)
526 WRITE (iw, *)
"Molecule kind:", i,
" contains", atom_in_kind,
" atoms"
527 IF (atom_in_kind <= 1) cycle
528 CALL find_boundary(atom_info%map_mol_typ, natom, first, last, i)
529 WRITE (iw, *)
"Boundary atoms:", first, last
530 cpassert(last - first + 1 == atom_in_kind)
531 max_mol_num = maxval(atom_info%map_mol_num(first:last))
532 WRITE (iw, *)
"Number of molecules of kind", i,
"is ::", max_mol_num
533 atom_in_mol = atom_in_kind/max_mol_num
534 WRITE (iw, *)
"Number of atoms per each molecule:", atom_in_mol
535 WRITE (iw, *)
"MAP_MOL_TYP::", atom_info%map_mol_typ(first:last)
536 WRITE (iw, *)
"MAP_MOL_NUM::", atom_info%map_mol_num(first:last)
537 WRITE (iw, *)
"MAP_MOL_RES::", atom_info%map_mol_res(first:last)
539 DO j = 1, max_mol_num
540 IF (count(atom_info%map_mol_num(first:last) == j) /= atom_in_mol)
THEN
541 WRITE (iw, *)
"molecule type:", i,
"molecule num:", j,
" has ", &
542 count(atom_info%map_mol_num(first:last) == j), &
543 " atoms instead of ", atom_in_mol,
" ."
544 CALL cp_abort(__location__, &
545 "Two molecules of the same kind have "// &
546 "been created with different numbers of atoms!")
552 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
553 CALL timestop(handle)
565 TYPE(mp_para_env_type),
POINTER :: para_env
566 TYPE(section_vals_type),
POINTER :: subsys_section
568 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_bond'
570 CHARACTER(LEN=2) :: upper_sym_1
571 INTEGER :: cbond, handle, handle2, i, iatm1, iatm2, iatom, ibond, idim, iw, j, jatom, k, &
572 n_bonds, n_heavy_bonds, n_hydr_bonds, n_rep, natom, npairs, output_unit
573 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bond_a, bond_b,
list, map_nb
574 INTEGER,
DIMENSION(:),
POINTER :: isolated_atoms, tmp_v
575 LOGICAL :: connectivity_ok, explicit, print_info
576 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: h_list
577 REAL(kind=
dp) :: bondparm_factor, cell_v(3), dr(3), &
578 ksign, my_maxrad, r2, r2_min, rbond, &
580 REAL(kind=
dp),
DIMENSION(1, 1) :: r_max, r_minsq
581 REAL(kind=
dp),
DIMENSION(:),
POINTER :: radius
582 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pbc_coord
585 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
586 TYPE(atomic_kind_type),
POINTER :: atomic_kind
588 TYPE(cp_logger_type),
POINTER :: logger
589 TYPE(fist_neighbor_type),
POINTER :: nonbonded
590 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
591 TYPE(section_vals_type),
POINTER :: bond_section, generate_section, &
594 NULLIFY (logger, particle_set, atomic_kind_set, nonbonded, bond_section, generate_section)
595 NULLIFY (isolated_atoms, tmp_v)
596 CALL timeset(routinen, handle)
600 extension=
".subsysLog")
602 ALLOCATE (isolated_atoms(0))
610 CALL reallocate(isolated_atoms, 1,
SIZE(isolated_atoms) +
SIZE(tmp_v))
611 isolated_atoms(
SIZE(isolated_atoms) -
SIZE(tmp_v) + 1:
SIZE(isolated_atoms)) = tmp_v
616 bondparm_factor =
topology%bondparm_factor
621 ALLOCATE (radius(natom))
622 ALLOCATE (
list(natom))
623 ALLOCATE (h_list(natom))
624 ALLOCATE (pbc_coord(3, natom))
626 CALL timeset(trim(routinen)//
"_1", handle2)
629 upper_sym_1 =
id2str(atom_info%id_element(iatom))
631 CALL get_ptable_info(symbol=upper_sym_1, covalent_radius=radius(iatom))
635 cpabort(
"Illegal bondparm_type")
637 IF (upper_sym_1 ==
"H ") h_list(iatom) = .true.
639 IF (any(isolated_atoms == iatom)) radius(iatom) = 0.0_dp
641 IF (iw > 0)
WRITE (iw,
'(T2,"GENERATE|",5X,A,T50,A5,T60,A,T69,F12.6)') &
642 "In topology_generate_bond :: iatom = ", upper_sym_1, &
643 "radius:", radius(iatom)
645 CALL timestop(handle2)
646 CALL timeset(trim(routinen)//
"_2", handle2)
649 ALLOCATE (atomic_kind_set(1))
652 my_maxrad = maxval(radius)*2.0_dp
653 atomic_kind => atomic_kind_set(1)
655 name=
"XXX", element_symbol=
"XXX", mass=0.0_dp, atom_list=
list)
658 IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT.
topology%molname_generated))
THEN
659 IF (output_unit > 0)
THEN
660 WRITE (output_unit,
'(T2,"GENERATE|",A)') &
661 " ERROR in connectivity generation!", &
662 " The THRESHOLD to select possible bonds is larger than the max. bondlength", &
663 " used to build the neighbors lists. Increase the BONDLENGTH_MAX parameter"
664 WRITE (output_unit,
'(T2,"GENERATE|",2(A,F11.6),A)') &
665 " Present THRESHOLD (", my_maxrad*bondparm_factor,
" )."// &
666 " Present BONDLENGTH_MAX (", r_max(1, 1),
" )"
671 particle_set(i)%atomic_kind => atomic_kind_set(1)
672 particle_set(i)%r(1) = atom_info%r(1, i)
673 particle_set(i)%r(2) = atom_info%r(2, i)
674 particle_set(i)%r(3) = atom_info%r(3, i)
675 pbc_coord(:, i) =
pbc(atom_info%r(:, i),
topology%cell)
679 CALL timestop(handle2)
680 CALL timeset(trim(routinen)//
"_3", handle2)
682 cell=
topology%cell, r_max=r_max, r_minsq=r_minsq, &
683 ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nonbonded, &
684 para_env=para_env, build_from_scratch=.true., geo_check=.true., &
685 mm_section=generate_section)
687 WRITE (iw,
'(T2,"GENERATE| Number of prescreened bonds (neighbors):",T71,I10)') &
688 nonbonded%neighbor_kind_pairs(1)%npairs
691 DO i = 1,
SIZE(nonbonded%neighbor_kind_pairs)
692 npairs = npairs + nonbonded%neighbor_kind_pairs(i)%npairs
694 ALLOCATE (bond_a(npairs))
695 ALLOCATE (bond_b(npairs))
696 ALLOCATE (map_nb(npairs))
698 DO j = 1,
SIZE(nonbonded%neighbor_kind_pairs)
699 DO i = 1, nonbonded%neighbor_kind_pairs(j)%npairs
701 bond_a(idim) = nonbonded%neighbor_kind_pairs(j)%list(1, i)
702 bond_b(idim) = nonbonded%neighbor_kind_pairs(j)%list(2, i)
706 CALL timestop(handle2)
707 CALL timeset(trim(routinen)//
"_4", handle2)
709 ALLOCATE (bond_list(natom))
711 ALLOCATE (bond_list(i)%array1(0))
712 ALLOCATE (bond_list(i)%array2(0))
714 CALL reorder_structure(bond_list, bond_a, bond_b, map_nb,
SIZE(bond_a))
721 CALL reallocate(conn_info%bond_a, 1, 1)
722 CALL reallocate(conn_info%bond_b, 1, 1)
723 connectivity_ok = .false.
727 atom_info%id_molname =
str2id(
s2s(
"TO_DEFINE_LATER"))
733 DO WHILE (.NOT. connectivity_ok)
737 IF (h_list(iatm1)) cycle
738 DO j = 1,
SIZE(bond_list(iatm1)%array1)
739 iatm2 = bond_list(iatm1)%array1(j)
740 IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) cycle
741 IF (h_list(iatm2) .OR. (iatm2 <= iatm1)) cycle
742 k = bond_list(iatm1)%array2(j)
743 ksign = sign(1.0_dp, real(k, kind=
dp))
745 cell_v = matmul(
topology%cell%hmat, &
746 REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, kind=
dp))
747 dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
748 r2 = dot_product(dr, dr)
749 IF (r2 <= r_minsq(1, 1))
THEN
750 CALL cp_abort(__location__, &
751 "bond distance between atoms less then the smallest distance provided "// &
752 "in input "//cp_to_string(tmp)//
" [bohr]")
755 IF ((any(isolated_atoms == iatm1)) .OR. (any(isolated_atoms == iatm2))) cycle
759 rbond = radius(iatm1) + radius(iatm2)
761 rbond = max(radius(iatm1), radius(iatm2))
764 rbond2 = rbond2*(bondparm_factor)**2
766 IF (r2 <= rbond2)
THEN
767 n_heavy_bonds = n_heavy_bonds + 1
768 CALL add_bonds_list(conn_info, iatm1, iatm2, n_heavy_bonds)
773 n_bonds = n_heavy_bonds
776 IF (output_unit > 0)
WRITE (output_unit, *)
778 IF (.NOT. h_list(iatm1)) cycle
779 r2_min = huge(0.0_dp)
782 DO j = 1,
SIZE(bond_list(iatm1)%array1)
783 iatm2 = bond_list(iatm1)%array1(j)
785 IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) cycle
786 IF (h_list(iatm2) .AND. (iatm2 <= iatm1)) cycle
788 IF ((any(isolated_atoms == iatm1)) .OR. (any(isolated_atoms == iatm2))) cycle
790 k = bond_list(iatm1)%array2(j)
791 ksign = sign(1.0_dp, real(k, kind=
dp))
793 cell_v = matmul(
topology%cell%hmat, &
794 REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, kind=
dp))
795 dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
796 r2 = dot_product(dr, dr)
797 IF (r2 <= r_minsq(1, 1))
THEN
798 CALL cp_abort(__location__, &
799 "bond distance between atoms less then the smallest distance provided "// &
800 "in input "//cp_to_string(tmp)//
" [bohr]")
802 IF (r2 <= r2_min)
THEN
807 IF (ibond == -1)
THEN
808 IF (output_unit > 0 .AND. print_info)
THEN
809 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,I10,A)') &
810 "WARNING:: No connections detected for Hydrogen - Atom Nr:", iatm1,
" !"
813 n_hydr_bonds = n_hydr_bonds + 1
814 n_bonds = n_bonds + 1
815 CALL add_bonds_list(conn_info, min(iatm1, ibond), max(iatm1, ibond), n_bonds)
818 IF (output_unit > 0)
THEN
819 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)') &
820 " Preliminary Number of Bonds generated:", n_bonds
824 CALL connectivity_external_control(section=bond_section, &
825 iarray1=conn_info%bond_a, &
826 iarray2=conn_info%bond_b, &
829 output_unit=output_unit)
831 CALL reallocate(conn_info%bond_a, 1, n_bonds)
832 CALL reallocate(conn_info%bond_b, 1, n_bonds)
836 IF (.NOT.
topology%reorder_atom)
THEN
838 IF (output_unit > 0)
WRITE (output_unit,
'(T2,"GENERATE|",A)') &
839 " Molecules names have been generated. Now reordering particle set in order to have ", &
840 " atoms belonging to the same molecule in a sequential order."
842 connectivity_ok = .true.
845 connectivity_ok = check_generate_mol(conn_info%bond_a, conn_info%bond_b, &
846 atom_info, bondparm_factor, output_unit)
848 IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT.
topology%molname_generated))
THEN
849 IF (output_unit > 0)
THEN
850 WRITE (output_unit,
'(T2,"GENERATE|",A)') &
851 " ERROR in connectivity generation!", &
852 " The THRESHOLD to select possible bonds is bigger than the MAX bondlength", &
853 " used to build the neighbors lists. Increase the BONDLENGTH_MAX patameter"
854 WRITE (output_unit,
'(T2,"GENERATE|",2(A,F11.6),A)') &
855 " Present THRESHOLD (", my_maxrad*bondparm_factor,
" )."// &
856 " Present BONDLENGTH_MAX (", r_max(1, 1),
" )"
861 IF (connectivity_ok .AND. (output_unit > 0))
THEN
862 WRITE (output_unit,
'(T2,"GENERATE|",A)') &
863 " Achieved consistency in connectivity generation."
866 CALL timestop(handle2)
867 CALL timeset(trim(routinen)//
"_6", handle2)
870 DEALLOCATE (bond_list(i)%array1)
871 DEALLOCATE (bond_list(i)%array2)
873 DEALLOCATE (bond_list)
874 DEALLOCATE (pbc_coord)
880 CALL timestop(handle2)
881 IF (output_unit > 0 .AND. n_bonds > 0)
THEN
882 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of Bonds generated:", &
885 CALL timeset(trim(routinen)//
"_7", handle2)
887 CALL reallocate(conn_info%c_bond_a, 1, 0)
888 CALL reallocate(conn_info%c_bond_b, 1, 0)
890 DO ibond = 1,
SIZE(conn_info%bond_a)
891 iatom = conn_info%bond_a(ibond)
892 jatom = conn_info%bond_b(ibond)
893 IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
894 (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
895 (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom)))
THEN
896 IF (iw > 0)
WRITE (iw, *)
" PARA_RES, bond between molecules atom ", &
899 CALL reallocate(conn_info%c_bond_a, 1, cbond)
900 CALL reallocate(conn_info%c_bond_b, 1, cbond)
901 conn_info%c_bond_a(cbond) = iatom
902 conn_info%c_bond_b(cbond) = jatom
904 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom))
THEN
905 cpabort(
"Bonds between different molecule types?")
910 CALL timestop(handle2)
911 DEALLOCATE (isolated_atoms)
912 CALL timestop(handle)
914 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
927 FUNCTION check_generate_mol(bond_a, bond_b, atom_info, bondparm_factor, output_unit) &
929 INTEGER,
DIMENSION(:),
POINTER :: bond_a, bond_b
931 REAL(kind=
dp),
INTENT(INOUT) :: bondparm_factor
932 INTEGER,
INTENT(IN) :: output_unit
935 CHARACTER(len=*),
PARAMETER :: routinen =
'check_generate_mol'
937 CHARACTER(LEN=10) :: ctmp1, ctmp2, ctmp3
938 INTEGER :: handle, i, idim, itype, j, mol_natom, &
940 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: mol_info_tmp
941 INTEGER,
DIMENSION(:),
POINTER :: mol_map, mol_map_o, wrk
942 INTEGER,
DIMENSION(:, :),
POINTER :: mol_info
943 LOGICAL,
DIMENSION(:),
POINTER :: icheck
946 CALL timeset(routinen, handle)
948 natom =
SIZE(atom_info%id_atmname)
949 ALLOCATE (bond_list(natom))
951 ALLOCATE (bond_list(i)%array1(0))
953 CALL reorder_structure(bond_list, bond_a, bond_b,
SIZE(bond_a))
954 ALLOCATE (mol_map(natom))
955 ALLOCATE (mol_map_o(natom))
956 ALLOCATE (wrk(natom))
959 mol_map(i) = atom_info%id_molname(i)
963 CALL sort(mol_map, natom, wrk)
977 ALLOCATE (mol_info_tmp(natom, 2))
982 mol_info_tmp(1, 1) = itype
984 IF (mol_map(i) /= itype)
THEN
987 mol_info_tmp(nsize, 1) = itype
988 mol_info_tmp(nsize - 1, 2) = idim
994 mol_info_tmp(nsize, 2) = idim
996 ALLOCATE (mol_info(nsize, 4))
997 mol_info(1:nsize, 1:2) = mol_info_tmp(1:nsize, 1:2)
998 DEALLOCATE (mol_info_tmp)
1005 ALLOCATE (icheck(natom))
1008 IF (icheck(i)) cycle
1009 itype = mol_map_o(i)
1012 DO j = 1,
SIZE(mol_info)
1013 IF (itype == mol_info(j, 1))
EXIT
1015 mol_info(j, 3) = mol_info(j, 3) + 1
1016 IF (mol_info(j, 4) == 0) mol_info(j, 4) = mol_natom
1017 IF (mol_info(j, 4) /= mol_natom)
THEN
1023 bondparm_factor = bondparm_factor*1.05_dp
1024 IF (output_unit < 0)
EXIT
1025 WRITE (output_unit,
'(/,T2,"GENERATE|",A)')
" WARNING in connectivity generation!"
1026 WRITE (output_unit,
'(T2,"GENERATE|",A)') &
1027 ' Two molecules/residues named ('//trim(
id2str(itype))//
') have different '// &
1032 WRITE (output_unit,
'(T2,"GENERATE|",A)')
' Molecule starting at position ('// &
1033 trim(ctmp1)//
') has Nr. <'//trim(ctmp2)// &
1034 '> of atoms.',
' while the other same molecules have Nr. <'// &
1035 trim(ctmp3)//
'> of atoms!'
1036 WRITE (output_unit,
'(T2,"GENERATE|",A)') &
1037 ' Increasing bondparm_factor by 1.05.. An error was found in the generated', &
1038 ' connectivity. Retry...'
1039 WRITE (output_unit,
'(T2,"GENERATE|",A,F11.6,A,/)') &
1040 " Present value of BONDPARM_FACTOR (", bondparm_factor,
" )."
1046 DEALLOCATE (mol_info)
1047 DEALLOCATE (mol_map)
1048 DEALLOCATE (mol_map_o)
1051 DEALLOCATE (bond_list(i)%array1)
1053 DEALLOCATE (bond_list)
1054 CALL timestop(handle)
1055 END FUNCTION check_generate_mol
1071 SUBROUTINE connectivity_external_control(section, Iarray1, Iarray2, Iarray3, Iarray4, nvar, &
1072 topology, output_unit, is_impr)
1073 TYPE(section_vals_type),
POINTER :: section
1074 INTEGER,
DIMENSION(:),
POINTER :: iarray1, iarray2
1075 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: iarray3, iarray4
1076 INTEGER,
INTENT(INOUT) :: nvar
1078 INTEGER,
INTENT(IN) :: output_unit
1079 LOGICAL,
INTENT(IN),
OPTIONAL :: is_impr
1081 CHARACTER(LEN=8) :: fmt
1082 INTEGER :: do_action, do_it, i, j, k, n_rep, &
1083 n_rep_val, natom, new_size, nsize
1084 INTEGER,
DIMENSION(:),
POINTER :: atlist, ilist1, ilist2, ilist3, ilist4
1085 LOGICAL :: explicit, ip3, ip4
1089 ip3 =
PRESENT(iarray3)
1090 ip4 =
PRESENT(iarray4)
1092 IF (ip3) nsize = nsize + 1
1093 IF (ip3 .AND. ip4) nsize = nsize + 1
1099 NULLIFY (ilist1, ilist2, ilist3, ilist4, atlist)
1100 ALLOCATE (ilist1(nvar))
1101 ALLOCATE (ilist2(nvar))
1102 ilist1 = iarray1(1:nvar)
1103 ilist2 = iarray2(1:nvar)
1107 ALLOCATE (ilist3(nvar))
1108 ilist3 = iarray3(1:nvar)
1110 ALLOCATE (ilist3(nvar))
1111 ALLOCATE (ilist4(nvar))
1112 ilist3 = iarray3(1:nvar)
1113 ilist4 = iarray4(1:nvar)
1118 CALL list_canonical_order(ilist1, ilist2, ilist3, ilist4, nsize, is_impr)
1128 cpassert(
SIZE(atlist) == nsize)
1130 CALL check_element_list(do_it, do_action, atlist, ilist1, ilist2, ilist3, ilist4, &
1132 IF (do_action ==
do_add)
THEN
1136 IF (output_unit > 0)
THEN
1137 WRITE (output_unit,
'(T2,"ADD|",1X,A,I6,'//trim(fmt)//
'(A,I6),A,T64,A,I6)') &
1139 atlist(1), (
",", atlist(k), k=2, nsize),
") added.",
" NEW size::", nvar
1141 IF (nvar >
SIZE(iarray1))
THEN
1142 new_size = int(5 + 1.2*nvar)
1143 CALL reallocate(iarray1, 1, new_size)
1144 CALL reallocate(iarray2, 1, new_size)
1147 CALL reallocate(iarray3, 1, new_size)
1149 CALL reallocate(iarray3, 1, new_size)
1150 CALL reallocate(iarray4, 1, new_size)
1154 iarray1(do_it + 1:nvar) = iarray1(do_it:nvar - 1)
1155 iarray2(do_it + 1:nvar) = iarray2(do_it:nvar - 1)
1156 iarray1(do_it) = ilist1(do_it)
1157 iarray2(do_it) = ilist2(do_it)
1160 iarray3(do_it + 1:nvar) = iarray3(do_it:nvar - 1)
1161 iarray3(do_it) = ilist3(do_it)
1163 iarray3(do_it + 1:nvar) = iarray3(do_it:nvar - 1)
1164 iarray4(do_it + 1:nvar) = iarray4(do_it:nvar - 1)
1165 iarray3(do_it) = ilist3(do_it)
1166 iarray4(do_it) = ilist4(do_it)
1169 IF (output_unit > 0)
THEN
1170 WRITE (output_unit,
'(T2,"ADD|",1X,A,I6,'//trim(fmt)//
'(A,I6),A,T80,A)') &
1172 atlist(1), (
",", atlist(k), k=2, nsize),
") already found.",
"X"
1179 IF (output_unit > 0)
THEN
1180 WRITE (output_unit,
'(T2,"RMV|",1X,A,I6,'//trim(fmt)//
'(A,I6),A,T64,A,I6)') &
1182 atlist(1), (
",", atlist(k), k=2, nsize),
") removed.",
" NEW size::", nvar
1184 iarray1(do_it:nvar) = iarray1(do_it + 1:nvar + 1)
1185 iarray2(do_it:nvar) = iarray2(do_it + 1:nvar + 1)
1186 iarray1(nvar + 1) = -huge(0)
1187 iarray2(nvar + 1) = -huge(0)
1190 iarray3(do_it:nvar) = iarray3(do_it + 1:nvar + 1)
1191 iarray3(nvar + 1) = -huge(0)
1193 iarray3(do_it:nvar) = iarray3(do_it + 1:nvar + 1)
1194 iarray4(do_it:nvar) = iarray4(do_it + 1:nvar + 1)
1195 iarray3(nvar + 1) = -huge(0)
1196 iarray4(nvar + 1) = -huge(0)
1199 IF (output_unit > 0)
THEN
1200 WRITE (output_unit,
'(T2,"RMV|",1X,A,I6,'//trim(fmt)//
'(A,I6),A,T80,A)') &
1202 atlist(1), (
",", atlist(k), k=2, nsize),
") not found.",
"X"
1223 END SUBROUTINE connectivity_external_control
1236 SUBROUTINE list_canonical_order(Ilist1, Ilist2, Ilist3, Ilist4, nsize, is_impr)
1237 INTEGER,
DIMENSION(:),
POINTER :: ilist1, ilist2
1238 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: ilist3, ilist4
1239 INTEGER,
INTENT(IN) :: nsize
1240 LOGICAL,
INTENT(IN),
OPTIONAL :: is_impr
1242 INTEGER :: i, ss(3), tmp1, tmp2, tmp3, tt(3)
1246 IF (
PRESENT(is_impr)) do_impr = is_impr
1249 DO i = 1,
SIZE(ilist1)
1252 ilist1(i) = min(tmp1, tmp2)
1253 ilist2(i) = max(tmp1, tmp2)
1256 DO i = 1,
SIZE(ilist1)
1259 ilist1(i) = min(tmp1, tmp2)
1260 ilist3(i) = max(tmp1, tmp2)
1263 DO i = 1,
SIZE(ilist1)
1264 IF (.NOT. do_impr)
THEN
1267 ilist1(i) = min(tmp1, tmp2)
1268 IF (ilist1(i) == tmp2)
THEN
1270 ilist3(i) = ilist2(i)
1273 ilist4(i) = max(tmp1, tmp2)
1278 CALL sort(tt, 3, ss)
1286 END SUBROUTINE list_canonical_order
1300 SUBROUTINE check_element_list(do_it, do_action, atlist, Ilist1, Ilist2, Ilist3, Ilist4, &
1302 INTEGER,
INTENT(OUT) :: do_it
1303 INTEGER,
INTENT(IN) :: do_action
1304 INTEGER,
DIMENSION(:),
POINTER :: atlist, ilist1, ilist2
1305 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: ilist3, ilist4
1306 LOGICAL,
INTENT(IN),
OPTIONAL :: is_impr
1308 INTEGER :: i, iend, istart, ndim, new_size, nsize, &
1309 ss(3), tmp1, tmp2, tmp3, tt(3)
1310 INTEGER,
DIMENSION(4) :: tmp
1311 LOGICAL :: do_impr, found
1314 IF (
PRESENT(is_impr)) do_impr = is_impr
1316 nsize =
SIZE(atlist)
1325 tmp(1) = min(tmp1, tmp2)
1326 tmp(2) = max(tmp1, tmp2)
1330 tmp(1) = min(tmp1, tmp2)
1331 tmp(3) = max(tmp1, tmp2)
1333 IF (.NOT. do_impr)
THEN
1336 tmp(1) = min(tmp1, tmp2)
1337 IF (tmp(1) == tmp2)
THEN
1342 tmp(4) = max(tmp1, tmp2)
1347 CALL sort(tt, 3, ss)
1355 IF (ilist1(istart) >= tmp(1))
EXIT
1358 IF (istart <= ndim)
THEN
1359 IF (ilist1(istart) > tmp(1) .AND. (istart /= 1)) istart = istart - 1
1361 DO iend = istart, ndim
1362 IF (ilist1(iend) /= tmp(1))
EXIT
1364 IF (iend == ndim + 1) iend = ndim
1369 IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2)))
EXIT
1370 IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2)))
THEN
1377 IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2)) .OR. (ilist3(i) > tmp(3)))
EXIT
1378 IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2)) .AND. (ilist3(i) == tmp(3)))
THEN
1385 IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2)) .OR. (ilist3(i) > tmp(3)) .OR. (ilist4(i) > tmp(4)))
EXIT
1386 IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2)) &
1387 .AND. (ilist3(i) == tmp(3)) .AND. (ilist4(i) == tmp(4)))
THEN
1393 SELECT CASE (do_action)
1408 CALL reallocate(ilist1, 1, new_size)
1409 CALL reallocate(ilist2, 1, new_size)
1410 ilist1(i + 1:new_size) = ilist1(i:ndim)
1411 ilist2(i + 1:new_size) = ilist2(i:ndim)
1416 CALL reallocate(ilist3, 1, new_size)
1417 ilist3(i + 1:new_size) = ilist3(i:ndim)
1420 CALL reallocate(ilist3, 1, new_size)
1421 CALL reallocate(ilist4, 1, new_size)
1422 ilist3(i + 1:new_size) = ilist3(i:ndim)
1423 ilist4(i + 1:new_size) = ilist4(i:ndim)
1433 ilist1(i:new_size) = ilist1(i + 1:ndim)
1434 ilist2(i:new_size) = ilist2(i + 1:ndim)
1435 CALL reallocate(ilist1, 1, new_size)
1436 CALL reallocate(ilist2, 1, new_size)
1439 ilist3(i:new_size) = ilist3(i + 1:ndim)
1440 CALL reallocate(ilist3, 1, new_size)
1442 ilist3(i:new_size) = ilist3(i + 1:ndim)
1443 ilist4(i:new_size) = ilist4(i + 1:ndim)
1444 CALL reallocate(ilist3, 1, new_size)
1445 CALL reallocate(ilist4, 1, new_size)
1454 END SUBROUTINE check_element_list
1464 SUBROUTINE add_bonds_list(conn_info, atm1, atm2, n_bonds)
1466 INTEGER,
INTENT(IN) :: atm1, atm2, n_bonds
1468 INTEGER :: new_size, old_size
1470 old_size =
SIZE(conn_info%bond_a)
1471 IF (n_bonds > old_size)
THEN
1472 new_size = int(5 + 1.2*old_size)
1473 CALL reallocate(conn_info%bond_a, 1, new_size)
1474 CALL reallocate(conn_info%bond_b, 1, new_size)
1476 conn_info%bond_a(n_bonds) = atm1
1477 conn_info%bond_b(n_bonds) = atm2
1478 END SUBROUTINE add_bonds_list
1488 TYPE(section_vals_type),
POINTER :: subsys_section
1490 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_bend'
1492 INTEGER :: handle, handle2, i, iw, natom, nbond, &
1493 nsize, ntheta, output_unit
1496 TYPE(cp_logger_type),
POINTER :: logger
1497 TYPE(section_vals_type),
POINTER :: bend_section
1502 extension=
".subsysLog")
1503 CALL timeset(routinen, handle)
1510 IF (
ASSOCIATED(conn_info%bond_a))
THEN
1511 nbond =
SIZE(conn_info%bond_a)
1513 CALL reallocate(conn_info%bond_a, 1, nbond)
1514 CALL reallocate(conn_info%bond_b, 1, nbond)
1516 IF (nbond /= 0)
THEN
1517 nsize = int(5 + 1.2*ntheta)
1518 CALL reallocate(conn_info%theta_a, 1, nsize)
1519 CALL reallocate(conn_info%theta_b, 1, nsize)
1520 CALL reallocate(conn_info%theta_c, 1, nsize)
1522 ALLOCATE (bond_list(natom))
1524 ALLOCATE (bond_list(i)%array1(0))
1526 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1528 CALL timeset(routinen//
"_1", handle2)
1529 CALL match_iterative_path(iarray1=bond_list, &
1530 iarray2=bond_list, &
1533 oarray1=conn_info%theta_a, &
1534 oarray2=conn_info%theta_b, &
1535 oarray3=conn_info%theta_c)
1536 CALL timestop(handle2)
1538 DEALLOCATE (bond_list(i)%array1)
1540 DEALLOCATE (bond_list)
1541 IF (output_unit > 0)
THEN
1542 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Preliminary Number of Bends generated:", &
1547 CALL connectivity_external_control(section=bend_section, &
1548 iarray1=conn_info%theta_a, &
1549 iarray2=conn_info%theta_b, &
1550 iarray3=conn_info%theta_c, &
1553 output_unit=output_unit)
1556 CALL reallocate(conn_info%theta_a, 1, ntheta)
1557 CALL reallocate(conn_info%theta_b, 1, ntheta)
1558 CALL reallocate(conn_info%theta_c, 1, ntheta)
1559 IF (output_unit > 0 .AND. ntheta > 0)
THEN
1560 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of Bends generated:", &
1563 CALL timestop(handle)
1565 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1585 RECURSIVE SUBROUTINE match_iterative_path(Iarray1, Iarray2, Iarray3, &
1586 max_levl, Oarray1, Oarray2, Oarray3, Oarray4, Ilist, it_levl, nvar)
1589 POINTER :: iarray2, iarray3
1590 INTEGER,
INTENT(IN) :: max_levl
1591 INTEGER,
DIMENSION(:),
POINTER :: oarray1, oarray2
1592 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: oarray3, oarray4
1593 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: ilist
1594 INTEGER,
INTENT(IN),
OPTIONAL :: it_levl
1595 INTEGER,
INTENT(INOUT) :: nvar
1597 INTEGER :: i, ind, j, my_levl, natom
1598 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: my_list
1602 check = max_levl >= 2 .AND. max_levl <= 4
1604 IF (.NOT.
PRESENT(ilist))
THEN
1605 SELECT CASE (max_levl)
1607 cpassert(.NOT.
PRESENT(iarray2))
1608 cpassert(.NOT.
PRESENT(iarray3))
1609 cpassert(.NOT.
PRESENT(oarray3))
1610 cpassert(.NOT.
PRESENT(oarray4))
1612 cpassert(
PRESENT(iarray2))
1613 cpassert(.NOT.
PRESENT(iarray3))
1614 cpassert(
PRESENT(oarray3))
1615 cpassert(.NOT.
PRESENT(oarray4))
1617 cpassert(
PRESENT(iarray2))
1618 cpassert(
PRESENT(iarray3))
1619 cpassert(
PRESENT(oarray3))
1620 cpassert(
PRESENT(oarray4))
1623 natom =
SIZE(iarray1)
1624 IF (.NOT.
PRESENT(ilist))
THEN
1626 ALLOCATE (my_list(max_levl))
1630 my_list(my_levl) = i
1631 CALL match_iterative_path(iarray1=iarray1, &
1634 it_levl=my_levl + 1, &
1635 max_levl=max_levl, &
1643 DEALLOCATE (my_list)
1645 SELECT CASE (it_levl)
1653 i = ilist(it_levl - 1)
1654 DO j = 1,
SIZE(iarray1(i)%array1)
1655 ind = wrk(i)%array1(j)
1656 IF (any(ilist == ind)) cycle
1657 IF (it_levl < max_levl)
THEN
1658 ilist(it_levl) = ind
1659 CALL match_iterative_path(iarray1=iarray1, &
1662 it_levl=it_levl + 1, &
1663 max_levl=max_levl, &
1671 ELSEIF (it_levl == max_levl)
THEN
1672 IF (ilist(1) > ind) cycle
1673 ilist(it_levl) = ind
1675 SELECT CASE (it_levl)
1677 IF (nvar >
SIZE(oarray1))
THEN
1678 CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1679 CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1681 oarray1(nvar) = ilist(1)
1682 oarray2(nvar) = ilist(2)
1684 IF (nvar >
SIZE(oarray1))
THEN
1685 CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1686 CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1687 CALL reallocate(oarray3, 1, int(5 + 1.2*nvar))
1689 oarray1(nvar) = ilist(1)
1690 oarray2(nvar) = ilist(2)
1691 oarray3(nvar) = ilist(3)
1693 IF (nvar >
SIZE(oarray1))
THEN
1694 CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1695 CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1696 CALL reallocate(oarray3, 1, int(5 + 1.2*nvar))
1697 CALL reallocate(oarray4, 1, int(5 + 1.2*nvar))
1699 oarray1(nvar) = ilist(1)
1700 oarray2(nvar) = ilist(2)
1701 oarray3(nvar) = ilist(3)
1702 oarray4(nvar) = ilist(4)
1714 END SUBROUTINE match_iterative_path
1725 TYPE(section_vals_type),
POINTER :: subsys_section
1727 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_ub'
1729 INTEGER :: handle, itheta, iw, ntheta, output_unit
1731 TYPE(cp_logger_type),
POINTER :: logger
1736 extension=
".subsysLog")
1738 CALL timeset(routinen, handle)
1740 ntheta =
SIZE(conn_info%theta_a)
1741 CALL reallocate(conn_info%ub_a, 1, ntheta)
1742 CALL reallocate(conn_info%ub_b, 1, ntheta)
1743 CALL reallocate(conn_info%ub_c, 1, ntheta)
1745 DO itheta = 1, ntheta
1746 conn_info%ub_a(itheta) = conn_info%theta_a(itheta)
1747 conn_info%ub_b(itheta) = conn_info%theta_b(itheta)
1748 conn_info%ub_c(itheta) = conn_info%theta_c(itheta)
1750 IF (output_unit > 0 .AND. ntheta > 0)
THEN
1751 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of UB generated:", &
1754 CALL timestop(handle)
1756 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1768 TYPE(section_vals_type),
POINTER :: subsys_section
1770 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_dihe'
1772 INTEGER :: handle, i, iw, natom, nbond, nphi, &
1776 TYPE(cp_logger_type),
POINTER :: logger
1777 TYPE(section_vals_type),
POINTER :: torsion_section
1782 extension=
".subsysLog")
1784 CALL timeset(routinen, handle)
1787 nbond =
SIZE(conn_info%bond_a)
1788 IF (nbond /= 0)
THEN
1789 nsize = int(5 + 1.2*nphi)
1790 CALL reallocate(conn_info%phi_a, 1, nsize)
1791 CALL reallocate(conn_info%phi_b, 1, nsize)
1792 CALL reallocate(conn_info%phi_c, 1, nsize)
1793 CALL reallocate(conn_info%phi_d, 1, nsize)
1796 ALLOCATE (bond_list(natom))
1798 ALLOCATE (bond_list(i)%array1(0))
1800 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1802 CALL match_iterative_path(iarray1=bond_list, &
1803 iarray2=bond_list, &
1804 iarray3=bond_list, &
1807 oarray1=conn_info%phi_a, &
1808 oarray2=conn_info%phi_b, &
1809 oarray3=conn_info%phi_c, &
1810 oarray4=conn_info%phi_d)
1812 DEALLOCATE (bond_list(i)%array1)
1814 DEALLOCATE (bond_list)
1815 IF (output_unit > 0)
THEN
1816 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Preliminary Number of Torsions generated:", &
1821 CALL connectivity_external_control(section=torsion_section, &
1822 iarray1=conn_info%phi_a, &
1823 iarray2=conn_info%phi_b, &
1824 iarray3=conn_info%phi_c, &
1825 iarray4=conn_info%phi_d, &
1828 output_unit=output_unit)
1831 CALL reallocate(conn_info%phi_a, 1, nphi)
1832 CALL reallocate(conn_info%phi_b, 1, nphi)
1833 CALL reallocate(conn_info%phi_c, 1, nphi)
1834 CALL reallocate(conn_info%phi_d, 1, nphi)
1835 IF (output_unit > 0 .AND. nphi > 0)
THEN
1836 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of Torsions generated:", &
1839 CALL timestop(handle)
1841 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1853 TYPE(section_vals_type),
POINTER :: subsys_section
1855 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_impr'
1857 CHARACTER(LEN=2) :: atm_symbol
1858 INTEGER :: handle, i, ind, iw, j, natom, nbond, &
1859 nimpr, nsize, output_unit
1860 LOGICAL :: accept_impr
1864 TYPE(cp_logger_type),
POINTER :: logger
1865 TYPE(section_vals_type),
POINTER :: impr_section
1870 extension=
".subsysLog")
1872 CALL timeset(routinen, handle)
1877 nbond =
SIZE(conn_info%bond_a)
1878 IF (nbond /= 0)
THEN
1879 nsize = int(5 + 1.2*nimpr)
1880 CALL reallocate(conn_info%impr_a, 1, nsize)
1881 CALL reallocate(conn_info%impr_b, 1, nsize)
1882 CALL reallocate(conn_info%impr_c, 1, nsize)
1883 CALL reallocate(conn_info%impr_d, 1, nsize)
1885 ALLOCATE (bond_list(natom))
1887 ALLOCATE (bond_list(i)%array1(0))
1889 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1892 IF (
SIZE(bond_list(i)%array1) == 3)
THEN
1895 accept_impr = .true.
1896 atm_symbol =
id2str(atom_info%id_element(i))
1898 IF (atm_symbol ==
"N ")
THEN
1899 accept_impr = .false.
1903 ind = bond_list(i)%array1(j)
1904 IF (
SIZE(bond_list(ind)%array1) == 3) accept_impr = .true.
1907 IF (.NOT. accept_impr) cycle
1909 IF (nimpr >
SIZE(conn_info%impr_a))
THEN
1910 nsize = int(5 + 1.2*nimpr)
1911 CALL reallocate(conn_info%impr_a, 1, nsize)
1912 CALL reallocate(conn_info%impr_b, 1, nsize)
1913 CALL reallocate(conn_info%impr_c, 1, nsize)
1914 CALL reallocate(conn_info%impr_d, 1, nsize)
1916 conn_info%impr_a(nimpr) = i
1917 conn_info%impr_b(nimpr) = bond_list(i)%array1(1)
1918 conn_info%impr_c(nimpr) = bond_list(i)%array1(2)
1919 conn_info%impr_d(nimpr) = bond_list(i)%array1(3)
1923 DEALLOCATE (bond_list(i)%array1)
1925 DEALLOCATE (bond_list)
1928 CALL connectivity_external_control(section=impr_section, &
1929 iarray1=conn_info%impr_a, &
1930 iarray2=conn_info%impr_b, &
1931 iarray3=conn_info%impr_c, &
1932 iarray4=conn_info%impr_d, &
1935 output_unit=output_unit, &
1939 CALL reallocate(conn_info%impr_a, 1, nimpr)
1940 CALL reallocate(conn_info%impr_b, 1, nimpr)
1941 CALL reallocate(conn_info%impr_c, 1, nimpr)
1942 CALL reallocate(conn_info%impr_d, 1, nimpr)
1943 IF (output_unit > 0 .AND. nimpr > 0)
THEN
1944 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of Impropers generated:", &
1947 CALL timestop(handle)
1949 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1960 TYPE(section_vals_type),
POINTER :: subsys_section
1962 CHARACTER(len=*),
PARAMETER :: routinen =
'topology_generate_onfo'
1964 INTEGER :: atom_a, atom_b, handle, i, ionfo, iw, &
1965 natom, nbond, nphi, ntheta, output_unit
1966 TYPE(
array1_list_type),
DIMENSION(:),
POINTER :: bond_list, phi_list, theta_list
1968 TYPE(cp_logger_type),
POINTER :: logger
1973 extension=
".subsysLog")
1975 CALL timeset(routinen, handle)
1981 ALLOCATE (bond_list(natom))
1983 ALLOCATE (bond_list(i)%array1(0))
1985 nbond =
SIZE(conn_info%bond_a)
1986 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1989 ALLOCATE (theta_list(natom))
1991 ALLOCATE (theta_list(i)%array1(0))
1993 ntheta =
SIZE(conn_info%theta_a)
1994 CALL reorder_structure(theta_list, conn_info%theta_a, conn_info%theta_c, ntheta)
1997 ALLOCATE (phi_list(natom))
1999 ALLOCATE (phi_list(i)%array1(0))
2001 nphi =
SIZE(conn_info%phi_a)
2002 CALL reorder_structure(phi_list, conn_info%phi_a, conn_info%phi_d, nphi)
2005 CALL reallocate(conn_info%onfo_a, 1, nphi)
2006 CALL reallocate(conn_info%onfo_b, 1, nphi)
2009 DO atom_a = 1, natom
2010 DO i = 1,
SIZE(phi_list(atom_a)%array1)
2011 atom_b = phi_list(atom_a)%array1(i)
2013 IF (atom_a > atom_b) cycle
2015 IF (any(atom_b == bond_list(atom_a)%array1)) cycle
2017 IF (any(atom_b == theta_list(atom_a)%array1)) cycle
2019 IF (any(atom_b == phi_list(atom_a)%array1(:i - 1))) cycle
2021 conn_info%onfo_a(ionfo) = atom_a
2022 conn_info%onfo_b(ionfo) = atom_b
2027 CALL reallocate(conn_info%onfo_a, 1, ionfo)
2028 CALL reallocate(conn_info%onfo_b, 1, ionfo)
2032 DEALLOCATE (bond_list(i)%array1)
2034 DEALLOCATE (bond_list)
2037 DEALLOCATE (theta_list(i)%array1)
2039 DEALLOCATE (theta_list)
2042 DEALLOCATE (phi_list(i)%array1)
2044 DEALLOCATE (phi_list)
2047 IF (output_unit > 0 .AND. ionfo > 0)
THEN
2048 WRITE (output_unit,
'(T2,"GENERATE|",1X,A,T71,I10)')
" Number of 1-4 interactions generated:", &
2051 CALL timestop(handle)
2053 "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Define the atomic kind types and their sub types.
subroutine, public set_atomic_kind(atomic_kind, element_symbol, name, mass, kind_number, natom, atom_list, fist_potential, shell, shell_active, damping)
Set the components of an atomic kind data set.
subroutine, public deallocate_atomic_kind_set(atomic_kind_set)
Destructor routine for a set of atomic kinds.
Handles all functions related to the CELL.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Define the neighbor list data types and the corresponding functionality.
subroutine, public fist_neighbor_deallocate(fist_neighbor)
...
Generate the atomic neighbor lists for FIST.
subroutine, public build_fist_neighbor_lists(atomic_kind_set, particle_set, local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, nonbonded, para_env, build_from_scratch, geo_check, mm_section, full_nl, exclusions)
...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Utility routines for the memory handling.
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public deallocate_particle_set(particle_set)
Deallocate a particle set.
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Periodic Table related data definitions.
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
generates a unique id number for a string (str2id) that can be used two compare two strings....
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Utilities for string manipulations.
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Collection of subroutine needed for topology related things.
subroutine, public topology_generate_impr(topology, subsys_section)
Using a list of bends, generate a list of impr.
subroutine, public topology_generate_onfo(topology, subsys_section)
Using a list of torsion, generate a list of onfo.
subroutine, public topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, id_molname)
Generates molnames: useful when the connectivity on file does not provide them.
subroutine, public topology_generate_bend(topology, subsys_section)
Using a list of bonds, generate a list of bends.
subroutine, public topology_generate_molecule(topology, qmmm, qmmm_env, subsys_section)
Use information from bond list to generate molecule. (ie clustering)
subroutine, public topology_generate_dihe(topology, subsys_section)
Generate a list of torsions from bonds.
subroutine, public topology_generate_ub(topology, subsys_section)
The list of Urey-Bradley is equal to the list of bends.
subroutine, public topology_generate_bond(topology, para_env, subsys_section)
Use info from periodic table and assumptions to generate bonds.
Collection of subroutine needed for topology related things.
subroutine, public find_molecule(atom_bond_list, mol_info, mol_name)
each atom will be assigned a molecule number based on bonded fragments The array mol_info should be i...
recursive subroutine, public give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
...
recursive subroutine, public reorder_list_array(Ilist1, Ilist2, Ilist3, Ilist4, nsize, ndim)
Order arrays of lists.
Control for reading in different topologies and coordinates.
All kind of helpful little routines.