PROGRESS  master
prg_system_mod.F90
Go to the documentation of this file.
1 
5 module prg_system_mod
6 
7  use bml
8  use iso_c_binding
9  use prg_openfiles_mod
10  use prg_ptable_mod
11  use prg_graph_mod
12  use prg_extras_mod
13 
14 #ifdef USE_NVTX
15  use prg_nvtx_mod
16 #endif
17 
18  implicit none
19 
20  private
21 
22  integer, parameter :: dp = kind(1.0d0)
23 
25  type, public :: estruct_type
26 
28  integer :: norbs
29 
31  !approaches
32  integer :: norbscore
33 
35  integer :: nel
36 
38  integer, allocatable :: hindex(:,:)
39 
41  type(bml_matrix_t) :: ham
42 
44  type(bml_matrix_t) :: hamaux
45 
47  type(bml_matrix_t) :: ham0
48 
50  type(bml_matrix_t) :: oham
51 
53  type(bml_matrix_t) :: ohamaux
54 
56  type(bml_matrix_t) :: over
57 
59  type(bml_matrix_t) :: rho
60 
62  type(bml_matrix_t) :: orho
63 
65  type(bml_matrix_t) :: zmat
66 
68  type(bml_matrix_t) :: evects
69 
71  type(bml_matrix_t) :: mat
72 
74  type(bml_matrix_t) :: mataux
75 
77  real(dp), allocatable :: coul_pot_r(:)
78 
80  real(dp), allocatable :: coul_pot_k(:)
81 
83  real(dp), allocatable :: skforce(:,:)
84 
86  real(dp), allocatable :: fpul(:,:)
87 
89  real(dp), allocatable :: fscoul(:,:)
90 
92  real(dp), allocatable :: aux(:,:)
93 
95  real(dp), allocatable :: evals(:)
96 
98  real(dp), allocatable :: dvals(:)
99 
101  real(dp), allocatable :: ker(:,:)
102 
104  integer :: norbsincore
105 
107  real(dp) :: eband
108 
109  end type estruct_type
110 
112  type, public :: system_type
113 
115  integer :: nats = 0
116 
123  character(2), allocatable :: symbol(:)
124 
126  integer, allocatable :: atomic_number(:)
127 
131  real(dp), allocatable :: coordinate(:,:)
132 
136  real(dp), allocatable :: velocity(:,:)
137 
141  real(dp), allocatable :: force(:,:)
142 
146  real(dp), allocatable :: net_charge(:)
147 
153  real(dp), allocatable :: mass(:)
154 
163  real(dp), allocatable :: lattice_vector(:,:)
164 
171  real(dp), allocatable :: recip_vector(:,:)
172 
175  real(dp) :: volr
176 
179  real(dp) :: volk
180 
186  integer :: nsp
187 
194  integer, allocatable :: spindex(:)
195 
202  character(2), allocatable :: splist(:)
203 
209  integer, allocatable :: spatnum(:)
210 
215  real(dp), allocatable :: spmass(:)
216 
218  real(dp), allocatable :: userdef(:)
219 
221  integer, allocatable :: resindex(:)
222 
224  character(3), allocatable :: resname(:)
225 
227  character(3), allocatable :: atomname(:)
228 
230  type(estruct_type) :: estr
231 
232  end type system_type
233 
234  public :: prg_parse_system, prg_get_nameandext, prg_make_random_system, prg_write_system, prg_write_trajectory
235  public :: prg_get_origin, prg_get_covgraph, prg_get_subsystem, prg_translateandfoldtobox, prg_molpartition, prg_get_partial_atomgraph
236  public :: prg_destroy_subsystems, prg_get_covgraph_h, prg_collect_graph_p, prg_merge_graph, prg_merge_graph_adj, prg_adj2bml, prg_graph2bml
237  public :: prg_graph2vector, prg_vector2graph, prg_sortadj, prg_get_recip_vects, prg_translatetogeomcandfoldtobox
238  public :: prg_write_trajectoryandproperty, prg_get_distancematrix, prg_parameters_to_vectors, prg_vectors_to_parameters
239  public :: prg_get_dihedral, prg_wraparound, prg_centeratbox, prg_replicate, prg_cleanuprepeatedatoms
240  public :: prg_replicate_system, prg_destroy_system, prg_destroy_estr
241  public :: prg_collect_extended_graph_p
242 
243 contains
244 
250  subroutine prg_get_nameandext(fullfilename,filename,ext)
251  implicit none
252  character(len=*), intent(in) :: fullfilename
253  character(50), intent(inout) :: filename
254  character(3), intent(inout) :: ext
255  character(1), allocatable :: tempc(:)
256  character(len=30) :: tempcflex
257  integer :: lenc
258 
259  lenc=len(adjustl(trim(fullfilename)))
260  if(.not.allocated(tempc))allocate(tempc(lenc))
261  tempcflex = adjustl(trim(fullfilename))
262  filename = adjustl(trim(tempcflex(1:lenc-4)))
263  ext = adjustl(trim(tempcflex(lenc-2:lenc+1)))
264 
265  end subroutine prg_get_nameandext
266 
267 
273  subroutine prg_parse_system(system,filename,extin)
274  implicit none
275  character(1) :: onechar
276  character(10) :: dummyc(10)
277  character(2) :: possiblenewsymbol, twochar
278  character(2), allocatable :: sptempsymbols(:)
279  character(3), optional, intent(in) :: extin
280  character(3) :: extension
281  character(30) :: dummy
282  character(50) :: io_name, nametmp
283  character(60) :: pdbformat
284  character(len=*) :: filename
285  integer :: dummyi(10), header_lines, i, io_unit
286  integer :: ios, itemp, j, lines_to_lattice
287  integer :: max_lines, nats, nsp, verbose
288  integer :: lines_to_masses, lines_to_atom
289  logical :: islattice
290  real(dp) :: abc_angles(2,3), dummyr(10), lbx1
291  real(dp) :: lbx2, lby1, lby2, lbz1
292  real(dp) :: lbz2, max_x, max_y, max_z
293  real(dp) :: min_x, min_y, min_z, scfactor
294  type(system_type), intent(out) :: system
295 
296  verbose = 0
297  islattice = .false.
298 
299  max_x = 0.0_dp; max_y = 0.0_dp; max_z = 0.0_dp;
300  min_x = 0.0_dp; min_y = 0.0_dp; min_z = 0.0_dp;
301 
302  if(allocated(system%symbol)) stop "ERROR: System already allocated"
303 
304  if(.not.present(extin))then
305  call prg_get_nameandext(filename,nametmp,extension)
306  else
307  extension = extin
308  nametmp = trim(adjustl(filename))
309  endif
310 
311  select case(extension)
312 
313  case("xyz")
314 
315  !! For xyz format see http://openbabel.org/wiki/XYZ_%28format%29
316  io_name=trim(adjustl(nametmp))//".xyz"
317  call prg_open_file_to_read(io_unit,io_name)
318  read(io_unit,*)nats
319  read(io_unit,*)
320  system%nats = nats
321  allocate(system%symbol(nats))
322  allocate(system%atomic_number(nats))
323  allocate(system%coordinate(3,nats))
324  allocate(system%mass(nats))
325  allocate(system%lattice_vector(3,3))
326 
327  system%lattice_vector = 0.0_dp
328 
329  do i=1,nats
330  read(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i)
331  system%atomic_number(i) = element_atomic_number(system%symbol(i))
332  system%mass(i) = element_mass(system%atomic_number(i))
333  max_x = max(max_x,system%coordinate(1,i))
334  min_x = min(min_x,system%coordinate(1,i))
335  max_y = max(max_y,system%coordinate(2,i))
336  min_y = min(min_y,system%coordinate(2,i))
337  max_z = max(max_z,system%coordinate(3,i))
338  min_z = min(min_z,system%coordinate(3,i))
339  ! write(*,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i)
340  enddo
341 
342  !The following is not part of an xyz format.
343  !VMD, babel, xmakemol and pymol can still read this file.
344  read(io_unit,*,iostat=ios)dummy
345  if(dummy.eq."#lattice")then
346  if(verbose.eq.1)write(*,*)"There is a lattice ..."
347  read(io_unit,*)system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
348  read(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
349  read(io_unit,*)system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
350  else
351  !We will create the lattice vectors if they are not pressent.
352  !It will add 10.0 Ang to each coordinate.
353  system%lattice_vector = 0.0_dp
354  system%lattice_vector(1,1) = max_x - min_x + 10.0_dp
355  system%lattice_vector(2,2) = max_y - min_y + 10.0_dp
356  system%lattice_vector(3,3) = max_z - min_z + 10.0_dp
357  endif
358 
359  close(io_unit)
360 
361  case("pdb")
362 
363  !! For PDB format see http://www.wwpdb.org/documentation/file-format
364  io_name=trim(adjustl(nametmp))//".pdb"
365  call prg_open_file_to_read(io_unit,io_name)
366  header_lines = 0
367  lines_to_lattice = 0
368  max_lines = 1000000
369 
370  !Counting header lines
371  do i=1,max_lines
372  read(io_unit,*)dummy
373  if(dummy.eq."ATOM".or.dummy.eq."HETATM") then
374  exit
375  else
376  header_lines = header_lines + 1
377  if(dummy.eq."CRYST1")then
378  lines_to_lattice = header_lines
379  islattice = .true.
380  endif
381  endif
382  enddo
383 
384  nats = 1
385  do i=1,max_lines
386  read(io_unit,*)dummy
387  if(dummy.eq."ATOM".or.dummy.eq."HETATM") then
388  nats = nats + 1
389  else
390  exit
391  endif
392  enddo
393  close(io_unit)
394 
395  call prg_open_file_to_read(io_unit,io_name)
396  do i=1,header_lines
397  if(i.eq.lines_to_lattice)then
398  read(io_unit,*)dummy,abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
399  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)
400  else
401  read(io_unit,*)dummy
402  endif
403  enddo
404 
405  system%nats = nats
406  allocate(system%symbol(nats))
407  allocate(system%atomic_number(nats))
408  allocate(system%coordinate(3,nats))
409  allocate(system%mass(nats))
410  allocate(system%lattice_vector(3,3))
411  allocate(system%resname(nats))
412  allocate(system%resindex(nats))
413  allocate(system%atomname(nats))
414 
415  system%lattice_vector = 0.0_dp
416 
417  pdbformat= '(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2)'
418 
419  do i=1,nats
420  read(io_unit,pdbformat)dummyc(1),dummyi(1), &
421  system%atomname(i),dummyc(3),system%resname(i),dummyc(4),system%resindex(i),dummyc(5),&
422  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
423  dummyr(1),dummyr(2),system%symbol(i),dummyc(10)
424 
425  ! In case there are no symbols in the last column:
426  if(dummyc(4).ne."".and.system%symbol(i).eq."")then
427  onechar=adjustl(trim(dummyc(4)))
428  if(onechar.ne."H".and. &
429  onechar.ne."B".and. &
430  onechar.ne."C".and. &
431  onechar.ne."N".and. &
432  onechar.ne."O".and. &
433  onechar.ne."F".and. &
434  onechar.ne."P".and. &
435  onechar.ne."S".and. &
436  onechar.ne."I")then
437  twochar=adjustl(trim(dummyc(4)))
438  system%symbol(i)=twochar
439  else
440  system%symbol(i)=onechar
441  endif
442  system%atomic_number(i) = element_atomic_number_upper(system%symbol(i))
443  system%symbol(i) = element_symbol(system%atomic_number(i)) !upper to lower char.
444  else
445  system%atomic_number(i) = element_atomic_number(system%symbol(i))
446  endif
447 
448  ! Getting the atomic mass
449  system%mass(i) = element_mass(system%atomic_number(i))
450 
451  ! Getting the system limits.
452  max_x = max(max_x,system%coordinate(1,i))
453  min_x = min(min_x,system%coordinate(1,i))
454  max_y = max(max_y,system%coordinate(2,i))
455  min_y = min(min_y,system%coordinate(2,i))
456  max_z = max(max_z,system%coordinate(3,i))
457  min_z = min(min_z,system%coordinate(3,i))
458 
459  enddo
460 
461  if(islattice)then
462  call prg_parameters_to_vectors(abc_angles,system%lattice_vector)
463  else
464  !We will create the lattice vectors if they are not pressent
465  !It will add 10.0 Ang to each coordinate.
466  system%lattice_vector = 0.0_dp
467  system%lattice_vector(1,1) = max_x - min_x + 10.0_dp
468  system%lattice_vector(2,2) = max_y - min_y + 10.0_dp
469  system%lattice_vector(3,3) = max_z - min_z + 10.0_dp
470  endif
471 
472  close(io_unit)
473 
474  case("ltt")
475 
476  !! For old inputblock.dat (old dat) format see the LATTE manual.
477  io_name=trim(adjustl(nametmp))//".ltt"
478  call prg_open_file_to_read(io_unit,io_name)
479  read(io_unit,*)dummy, nats
480  read(io_unit,*)scfactor
481  system%nats = nats
482  allocate(system%symbol(nats))
483  allocate(system%atomic_number(nats))
484  allocate(system%coordinate(3,nats))
485  allocate(system%mass(nats))
486  allocate(system%lattice_vector(3,3))
487 
488  !!The gets the llattice vectors from the box boundaries.
489  read(io_unit,*)lbx1,lbx2,lby1,lby2,lbz1,lbz2
490  system%lattice_vector = 0.0_dp
491  system%lattice_vector(1,1) = (lbx2-lbx1)*scfactor
492  system%lattice_vector(2,2) = (lby2-lby1)*scfactor
493  system%lattice_vector(3,3) = (lbz2-lbz1)*scfactor
494 
495  do i=1,nats
496  read(io_unit,*)system%symbol(i),system%coordinate(1,i)&
497  ,system%coordinate(2,i),system%coordinate(3,i)
498  system%coordinate(1,i)= system%coordinate(1,i)*scfactor
499  system%coordinate(2,i)= system%coordinate(2,i)*scfactor
500  system%coordinate(3,i)= system%coordinate(3,i)*scfactor
501  system%atomic_number(i) = element_atomic_number(system%symbol(i))
502  system%mass(i) = element_mass(system%atomic_number(i))
503  enddo
504 
505  close(io_unit)
506 
507  case("dat")
508 
509  !! For new inputblock.dat (dat) format see the LATTE manual.
510  io_name=trim(adjustl(nametmp))//".dat"
511  call prg_open_file_to_read(io_unit,io_name)
512  read(io_unit,*)nats
513  system%nats = nats
514  allocate(system%symbol(nats))
515  allocate(system%atomic_number(nats))
516  allocate(system%coordinate(3,nats))
517  allocate(system%mass(nats))
518  allocate(system%lattice_vector(3,3))
519 
520  !!The gets the llattice vectors from the box boundaries.
521  system%lattice_vector = 0.0_dp
522  read(io_unit,*)system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
523  read(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
524  read(io_unit,*)system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
525 
526  do i=1,nats
527  read(io_unit,*)system%symbol(i),system%coordinate(1,i)&
528  ,system%coordinate(2,i),system%coordinate(3,i)
529  system%coordinate(1,i)= system%coordinate(1,i)
530  system%coordinate(2,i)= system%coordinate(2,i)
531  system%coordinate(3,i)= system%coordinate(3,i)
532  system%atomic_number(i) = element_atomic_number(system%symbol(i))
533  system%mass(i) = element_mass(system%atomic_number(i))
534  enddo
535 
536  close(io_unit)
537 
538  case("gen")
539 
540  stop "gen format is not implemented for reading yet"
541 
542  case("lmp")
543 
544  !! For Lammps data.* input file.
545  !! For more information see: http://lammps.sandia.gov/doc/2001/data_format.html
546  io_name=adjustl(trim(nametmp))//".lmp"
547  call prg_open_file(io_unit,io_name)
548  read(io_unit,*)dummyc(1)
549  read(io_unit,*)system%nats
550 
551  do i=1,100
552  read(io_unit,*)dummyc(1),dummyc(2)!,dummyc(3)
553 
554  if(adjustl(trim(dummyc(2))) == "atom")then
555  lines_to_atom = i + 2
556  endif
557  if(adjustl(trim(dummyc(1))) == "Masses")then
558  lines_to_masses = i + 2
559  exit
560  endif
561  enddo
562  close(io_unit)
563 
564  call prg_open_file(io_unit,io_name)
565 
566  do i=1,lines_to_atom-1
567  read(io_unit,*)dummyc(1)
568  enddo
569 
570  read(io_unit,*)system%nsp
571 
572  do i=1,lines_to_masses - lines_to_atom - 4
573  read(io_unit,*)dummyc(1)
574  enddo
575 
576  read(io_unit,*)min_x, max_x
577  read(io_unit,*)min_y, max_y
578  read(io_unit,*)min_z, max_z
579 
580  allocate(system%lattice_vector(3,3))
581  system%lattice_vector = 0.0_dp
582  system%lattice_vector(1,1) = max_x - min_x
583  system%lattice_vector(2,2) = max_y - min_y
584  system%lattice_vector(3,3) = max_z - min_z
585 
586  read(io_unit,*)dummyc(1)
587 
588  allocate(system%spmass(system%nsp))
589  allocate(system%splist(system%nsp))
590 
591  do i=1,system%nsp
592  read(io_unit,*)dummyi(1),system%spmass(i)
593  enddo
594 
595  do i=1,system%nsp
596  do j=1,size(element_mass,dim=1)
597  if(abs(system%spmass(i) - element_mass(j)) < 1.0) then
598  system%splist(i) = element_symbol(j)
599  exit
600  endif
601  enddo
602  enddo
603 
604  do i=1,100000
605  read(io_unit,*)dummyc(1)
606  if(dummyc(1) == "Atoms")exit
607  enddo
608 
609  allocate(system%spindex(system%nats))
610  allocate(system%coordinate(3,system%nats))
611  allocate(system%symbol(system%nats))
612 
613  do i=1,system%nats
614  read(io_unit,*)dummyi(1),dummyi(2),system%spindex(i),dummyr(1),system%coordinate(1,i)&
615  ,system%coordinate(2,i),system%coordinate(3,i)
616  write(*,*)dummyi(1),dummyi(2),system%spindex(i),dummyr(1),system%coordinate(1,i)&
617  ,system%coordinate(2,i),system%coordinate(3,i)
618  system%symbol(i) = system%splist(system%spindex(i))
619  enddo
620 
621  close(io_unit)
622 
623  case default
624 
625  stop "The file extension is not valid. Only xyz, lmp, dat and pdb formats are implemented"
626 
627  end select
628 
629  !!Computing species and species index
630  !!Assigning species index to the system elements.
631  nsp=1
632  allocate(sptempsymbols(150)) !Temporary vector to accumulate the species.
633  sptempsymbols(1)=system%Symbol(1)
634 
635  if(.not.allocated(system%spindex)) allocate(system%spindex(nats))
636 
637  do i=1,system%nats
638  itemp = 0
639  do j=1,nsp !If there is a new chemical symbol
640  if(trim(system%symbol(i)).ne.trim(sptempsymbols(j)))then
641  itemp = itemp + 1
642  possiblenewsymbol = system%symbol(i)
643  endif
644  enddo
645  if(itemp == nsp)then
646  nsp = nsp+1
647  sptempsymbols(nsp)=possiblenewsymbol
648  endif
649  enddo
650 
651  if(allocated(system%splist))deallocate(system%splist)
652  if(allocated(system%spmass))deallocate(system%spmass)
653 
654  allocate(system%splist(nsp))
655  allocate(system%spmass(nsp))
656  allocate(system%spatnum(nsp))
657 
658  do i=1,nsp
659  system%splist(i) = sptempsymbols(i)
660  system%spatnum(i) = element_atomic_number(system%splist(i))
661  system%spmass(i) = element_mass(system%spatnum(i))
662  enddo
663 
664  deallocate(sptempsymbols)
665  system%nsp = nsp
666 
669  do i = 1,system%nats
670  do j=1,nsp
671  if(trim(system%symbol(i)).eq.trim(system%splist(j)))then
672  system%spindex(i)=j
673  endif
674  enddo
675  enddo
676 
677  end subroutine prg_parse_system
678 
679 
683  subroutine prg_destroy_system(sy)
684  implicit none
685  type(system_type), intent(inout) :: sy
686 
687  call prg_destroy_estr(sy%estr)
688  if(allocated(sy%symbol)) deallocate(sy%symbol)
689  if(allocated(sy%atomic_number)) deallocate(sy%atomic_number)
690  if(allocated(sy%coordinate)) deallocate(sy%coordinate)
691  if(allocated(sy%velocity)) deallocate(sy%velocity)
692  if(allocated(sy%force)) deallocate(sy%force)
693  if(allocated(sy%net_charge)) deallocate(sy%net_charge)
694  if(allocated(sy%mass)) deallocate(sy%mass)
695  if(allocated(sy%lattice_vector)) deallocate(sy%lattice_vector)
696  if(allocated(sy%recip_vector)) deallocate(sy%recip_vector)
697  if(allocated(sy%spindex)) deallocate(sy%spindex)
698  if(allocated(sy%splist)) deallocate(sy%splist)
699  if(allocated(sy%spatnum)) deallocate(sy%spatnum)
700  if(allocated(sy%spmass)) deallocate(sy%spmass)
701  if(allocated(sy%userdef)) deallocate(sy%userdef)
702  if(allocated(sy%resindex)) deallocate(sy%resindex)
703  if(allocated(sy%resname)) deallocate(sy%resname)
704  if(allocated(sy%atomname)) deallocate(sy%atomname)
705 
706  end subroutine prg_destroy_system
707 
711  subroutine prg_destroy_estr(estr)
712  implicit none
713  type(estruct_type), intent(inout) :: estr
714 
715  if(allocated(estr%hindex))deallocate(estr%hindex)
716  if(bml_allocated(estr%ham))call bml_deallocate(estr%ham)
717  if(bml_allocated(estr%hamAux)) call bml_deallocate(estr%hamAux)
718  if(bml_allocated(estr%ham0)) call bml_deallocate(estr%ham0)
719  if(bml_allocated(estr%oham)) call bml_deallocate(estr%oham)
720  if(bml_allocated(estr%ohamAux)) call bml_deallocate(estr%ohamAux)
721  if(bml_allocated(estr%over)) call bml_deallocate(estr%over)
722  if(bml_allocated(estr%rho)) call bml_deallocate(estr%rho)
723  if(bml_allocated(estr%evects)) call bml_deallocate(estr%evects)
724  if(bml_allocated(estr%matAux)) call bml_deallocate(estr%matAux)
725  if(bml_allocated(estr%zmat)) call bml_deallocate(estr%zmat)
726  if(allocated(estr%coul_pot_r)) deallocate(estr%coul_pot_r)
727  if(allocated(estr%coul_pot_k)) deallocate(estr%coul_pot_k)
728  if(allocated(estr%skforce)) deallocate(estr%skforce)
729  if(allocated(estr%fpul)) deallocate(estr%fpul)
730  if(allocated(estr%fscoul)) deallocate(estr%fscoul)
731  if(allocated(estr%aux)) deallocate(estr%aux)
732  if(allocated(estr%evals)) deallocate(estr%evals)
733  if(allocated(estr%dvals)) deallocate(estr%dvals)
734  if(allocated(estr%ker)) deallocate(estr%ker)
735 
736  end subroutine prg_destroy_estr
737 
738 
744  subroutine prg_write_system(system,filename,extin)
745  implicit none
746  character(len=*) :: filename
747  character(10) :: dummyc(10)
748  character(34) :: xyzformat
749  character(3), optional, intent(in) :: extin
750  character(3) :: extension
751  character(50) :: io_name, nametmp
752  character(60) :: pdbformat
753  character(100) :: io_message
754  integer :: dummyi(10), i, io_unit, nats
755  real(dp) :: abc_angles(2,3), dummyr(10), origin(3)
756  type(system_type), intent(inout) :: system
757 
758  nats = system%nats
759 
760  if(.not.present(extin))then
761  call prg_get_nameandext(filename,nametmp,extension)
762  else
763  extension = extin
764  nametmp = trim(filename)
765  endif
766 
767  select case(extension)
768 
769  case("xyz")
770 
771  io_name=trim(nametmp)//".xyz"
772  call prg_open_file(io_unit,io_name)
773  write(io_unit,*)nats
774  io_message = trim(adjustl(io_name))//" Generated by the PROGRESS library"
775  write(io_unit,*)trim(adjustl(io_message))
776  xyzformat = '(A2,A2,1F10.5,A2,1F10.5,A2,1F10.5)'
777  do i=1,nats
778  write(io_unit,xyzformat)system%symbol(i)," ",system%coordinate(1,i)," ",system%coordinate(2,i)," ",system%coordinate(3,i)
779  enddo
780 
781  !The following is not part of an xyz format but
782  !VMD, babel, xmakemol and pymol can still read this file.
783  write(io_unit,*)"#lattice vectors"
784  write(io_unit,"(3F20.5)")system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
785  write(io_unit,"(3F20.5)")system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
786  write(io_unit,"(3F20.5)")system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
787 
788  close(io_unit)
789 
790  case("pdb")
791 
792  dummyi = 0
793  dummyc = ""
794  dummyc(1) = "ATOM"
795  dummyr = 0.0_dp
796 
797  io_name=trim(nametmp)//".pdb"
798  call prg_open_file(io_unit,io_name)
799 
800  write(io_unit,'(A6,1X,A50)')"REMARK","Generated by PROGRESS library"
801  write(io_unit,'(A5,1X,A20)')"TITLE",io_name
802 
803  if(allocated(system%lattice_vector))then
804  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
805  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
806  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
807  endif
808 
809  write(io_unit,'(A5,A20)')"MODEL"," 1"
810 
811  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2 )'
812 
813  if(.not.allocated(system%resindex))then
814  allocate(system%resindex(nats))
815  system%resindex = 1
816  endif
817 
818  if(.not.allocated(system%resname))then
819  allocate(system%resname(nats))
820  system%resname = "MOL"
821  endif
822 
823  if(.not.allocated(system%atomname))then
824  allocate(system%atomname(nats))
825  system%atomname = system%symbol
826  endif
827 
828  if(allocated(system%net_charge))then
829  do i=1,nats
830  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
831  system%atomname(i),dummyc(5),system%resname(i),dummyc(7),system%resindex(i),dummyc(8),&
832  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
833  dummyr(1),system%net_charge(i),system%symbol(i),dummyc(10)
834  enddo
835  else
836  do i=1,nats
837  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
838  system%atomname(i),dummyc(5),system%resname(i),dummyc(7),system%resindex(i),dummyc(8),&
839  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
840  dummyr(1),dummyr(2),system%symbol(i),dummyc(10)
841  enddo
842  endif
843 
844  write(io_unit,'(A3)')"TER"
845  write(io_unit,'(A3)')"ENDMDL"
846  close(io_unit)
847 
848 
849  case("ltt")
850 
851  !! For old inputblock.dat (old dat) format see the LATTE manual.
852  io_name=trim(nametmp)//".ltt"
853  call prg_open_file(io_unit,io_name)
854  write(io_unit,*)"NATS= ", system%nats
855 
856  !Here we will always write the system in 1:1 scale
857  write(io_unit,*)"1.0"," Generated by the PROGRESS library"
858 
859  !! This gets the lattice vectors from the box boundaries.
860  write(io_unit,"(6F10.5)")0.0,system%lattice_vector(1,1),0.0,system%lattice_vector(2,2)&
861  ,0.0,system%lattice_vector(3,3)
862 
863  do i=1,nats
864  write(io_unit,"(A2,3F10.5)")system%symbol(i),system%coordinate(1,i)&
865  ,system%coordinate(2,i),system%coordinate(3,i)
866  enddo
867 
868  close(io_unit)
869 
870  case("dat")
871 
872  !! For new inputblock.dat (nat) format see the LATTE manual.
873  io_name=trim(nametmp)//".dat"
874  call prg_open_file(io_unit,io_name)
875  write(io_unit,*)system%nats
876 
877  !! This gets the lattice vectors from the box boundaries.
878  write(io_unit,"(3F10.5)")system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
879  write(io_unit,"(3F10.5)")system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
880  write(io_unit,"(3F10.5)")system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
881 
882  do i=1,nats
883  write(io_unit,"(A2,3F10.5)")system%symbol(i),system%coordinate(1,i)&
884  ,system%coordinate(2,i),system%coordinate(3,i)
885  enddo
886 
887  close(io_unit)
888 
889  case("gen")
890 
891  !! For gen format see DFTB+ manual.
892  io_name=trim(nametmp)//".gen"
893  call prg_open_file(io_unit,io_name)
894  write(io_unit,*)system%nats,"S"
895  write(io_unit,'(103A3)')(system%splist(i),i=1,system%nsp)
896  do i=1,nats
897  write(io_unit,"(I10,I10,3F10.5)")i,system%spindex(i),system%coordinate(1,i)&
898  ,system%coordinate(2,i),system%coordinate(3,i)
899  enddo
900  write(io_unit,*)"0.0 0.0 0.0"
901  write(io_unit,'(3F10.5)')system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
902  write(io_unit,'(3F10.5)')system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
903  write(io_unit,'(3F10.5)')system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
904 
905  close(io_unit)
906 
907  case("lmp")
908 
909  !! For Lammps data.* input file.
910  !! For more information see: http://lammps.sandia.gov/doc/2001/data_format.html
911  io_name= adjustl(trim(nametmp))//".lmp"
912  call prg_open_file(io_unit,io_name)
913  write(io_unit,*)"LAMMPS Description"
914  write(io_unit,*)""
915  write(io_unit,*)system%nats,"atoms"
916  write(io_unit,*)""
917  write(io_unit,*)system%nsp,"atom types"
918  write(io_unit,*)""
919  write(io_unit,*)0.0_dp,system%lattice_vector(1,1),"xlo xhi"
920  write(io_unit,*)0.0_dp,system%lattice_vector(2,2),"ylo yhi"
921  write(io_unit,*)0.0_dp,system%lattice_vector(3,3),"zlo zhi"
922  write(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(3,1), &
923  &system%lattice_vector(3,2), "xy xz yz"
924 
925  write(io_unit,*)""
926  write(io_unit,*)"Masses"
927  write(io_unit,*)""
928  do i=1,system%nsp
929  write(io_unit,*)" ",i,system%spmass(i)
930  enddo
931  write(io_unit,*)""
932  write(io_unit,*)"Atoms"
933  write(io_unit,*)""
934  do i=1,nats
935  write(io_unit,'(3I5,A6,3F10.5)')i,1,system%spindex(i)," 0.0",system%coordinate(1,i)&
936  ,system%coordinate(2,i),system%coordinate(3,i)
937  enddo
938 
939  close(io_unit)
940 
941  case default
942 
943  stop "The file extension is not valid. Only xyz, dat, lmp, pdb or gen formats are implemented"
944 
945  end select
946 
947  end subroutine prg_write_system
948 
956  subroutine prg_write_trajectory(system,iter,each,prg_deltat,filename,extension)
957  implicit none
958  character(len=*) :: filename
959  character(10) :: dummyc(10)
960  character(11) :: xyzformat
961  character(20) :: io_name
962  character(3) :: extension
963  character(60) :: pdbformat
964  integer :: dummyi(10), i, io_unit, nats
965  integer, intent(in) :: iter, each
966  integer, allocatable :: resindex(:)
967  real(dp) :: abc_angles(2,3), dummyr(10)
968  real(dp), intent(in) :: prg_deltat
969  type(system_type), intent(in) :: system
970 
971  if(mod(iter,each).eq.0.or.iter.eq.0.or.iter.eq.1)then
972  nats = system%nats
973 
974  select case(extension)
975 
976  case("xyz")
977 
978  io_name=trim(filename)//".xyz"
979  io_unit=get_file_unit(-1)
980 
981  if(iter.eq.0.or.iter.eq.1)then
982  open(unit=io_unit,file=io_name,status='unknown')
983  else
984  open(unit=io_unit,file=io_name,position = 'append',status='old')
985  endif
986 
987  write(io_unit,*)system%nats
988  write(io_unit,*)"frame", iter
989  if(allocated(system%net_charge))then
990  do i=1,nats
991  write(io_unit,"(A3,3X,F20.15,3X,F20.15,3X,F20.15,3X,F20.15)")system%symbol(i),system%coordinate(1,i),&
992  system%coordinate(2,i),system%coordinate(3,i),system%net_charge(i)
993  enddo
994  else
995  do i=1,nats
996  write(io_unit,"(A3,3X,F20.15,3X,F20.15,3X,F20.15)")system%symbol(i),system%coordinate(1,i),&
997  system%coordinate(2,i),system%coordinate(3,i)
998  enddo
999  endif
1000 
1001  close(io_unit)
1002 
1003  case("pdb")
1004 
1005  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
1006 
1007  dummyi = 0
1008  dummyc = ""
1009  dummyc(1) = "ATOM"
1010  dummyr = 0.0_dp
1011 
1012  io_name=trim(filename)//".pdb"
1013  io_unit=get_file_unit(-1)
1014 
1015  if(iter.eq.0.or.iter.eq.1)then
1016  open(unit=io_unit,file=io_name,status='unknown')
1017  else
1018  open(unit=io_unit,file=io_name,position = 'append',status='old')
1019  endif
1020 
1021  write(io_unit,'(A17,A4,F10.5)')"REMARK Trajectory"," t=",iter*prg_deltat
1022  write(io_unit,'(A31)')"REMARK THIS IS A SIMULATION BOX"
1023  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
1024  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
1025  write(io_unit,'(A5,I6)')"MODEL",iter
1026 
1027  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,11X,A2,A2 )'
1028 
1029  if(.not.allocated(system%resindex))then
1030  allocate(resindex(nats))
1031  resindex = 1
1032  else
1033  resindex = system%resindex
1034  endif
1035 
1036 
1037  if(allocated(system%net_charge))then
1038  do i=1,nats
1039  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1040  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1041  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1042  dummyr(1),system%net_charge(i),system%symbol(i),dummyc(10)
1043  enddo
1044  else
1045  do i=1,nats
1046  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1047  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1048  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1049  dummyr(1),dummyr(2),system%symbol(i),dummyc(10)
1050  enddo
1051  endif
1052 
1053  write(io_unit,'(A3)')"TER"
1054  write(io_unit,'(A6)')"ENDMDL"
1055  close(io_unit)
1056 
1057  case("dat")
1058 
1059  write(*,*)"write traj not implemented for dat "
1060 
1061  case default
1062 
1063  stop "The file extension is not valid. Only pdb and xyz formats are implemented"
1064 
1065  end select
1066 
1067  endif
1068 
1069  end subroutine prg_write_trajectory
1070 
1080  subroutine prg_write_trajectoryandproperty(system,iter,each,prg_deltat,scalarprop,filename,extension)
1081  implicit none
1082  character(len=*) :: filename
1083  character(10) :: dummyc(10)
1084  character(11) :: xyzformat
1085  character(20) :: io_name
1086  character(3) :: extension
1087  character(60) :: pdbformat
1088  integer :: dummyi(10), i, io_unit, nats
1089  integer, intent(in) :: iter, each
1090  integer, allocatable :: resindex(:)
1091  real(dp) :: abc_angles(2,3), dummyr(10)
1092  real(dp) :: maxprop, minprop, realtmp
1093  real(dp), intent(in) :: prg_deltat, scalarprop(:)
1094  type(system_type), intent(in) :: system
1095 
1096 
1097  ! maxprop = maxval(scalarprop(:))
1098  ! minprop = minval(scalarprop(:))
1099 
1100  if(mod(iter,each).eq.0.or.iter.eq.0.or.iter.eq.1)then
1101  nats = system%nats
1102 
1103  select case(extension)
1104 
1105  case("pdb")
1106 
1107  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
1108 
1109  dummyi = 0
1110  dummyc = ""
1111  dummyc(1) = "ATOM"
1112  dummyr = 0.0_dp
1113 
1114  io_name=trim(filename)//".pdb"
1115  io_unit=get_file_unit(-1)
1116 
1117  if(iter.eq.0.or.iter.eq.1)then
1118  open(unit=io_unit,file=io_name,status='unknown')
1119  else
1120  open(unit=io_unit,file=io_name,position = 'append',status='old')
1121  endif
1122 
1123  write(io_unit,'(A8,A4,F10.5)')"Trajectory"," t=",iter*prg_deltat
1124  write(io_unit,'(A24)')"THIS IS A SIMULATION BOX"
1125  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
1126  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
1127  write(io_unit,'(A5,I6)')"MODEL",iter
1128 
1129  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2 )'
1130 
1131  if(.not.allocated(system%resindex))then
1132  allocate(resindex(nats))
1133  resindex = 1
1134  else
1135  resindex = system%resindex
1136  endif
1137 
1138  do i=1,nats
1139  ! realtmp = (scalarprop(i) - minprop)/(maxprop-minprop) - 1.0_dp !Normalization
1140  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1141  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1142  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1143  dummyr(1),scalarprop(i),system%symbol(i),dummyc(10)
1144  enddo
1145 
1146  write(io_unit,'(A3)')"TER"
1147  write(io_unit,'(A3)')"ENDMDL"
1148  close(io_unit)
1149 
1150  case("xyz")
1151 
1152  io_name=trim(filename)//".xyz"
1153  io_unit=get_file_unit(-1)
1154 
1155  if(iter.eq.0.or.iter.eq.1)then
1156  open(unit=io_unit,file=io_name,status='unknown')
1157  else
1158  open(unit=io_unit,file=io_name,position = 'append',status='old')
1159  endif
1160 
1161  write(io_unit,*)system%nats
1162  write(io_unit,*)"frame", iter
1163  do i=1,nats
1164  ! realtmp = (scalarprop(i) - minprop)/(maxprop-minprop) - 1.0_dp !Normalization
1165  write(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1166  scalarprop(i)
1167  ! write(*,*)scalarprop(i)
1168  enddo
1169 
1170  close(io_unit)
1171 
1172  case default
1173 
1174  stop "The file extension is not valid. Only pdb format is implemented"
1175 
1176  end select
1177 
1178  endif
1179 
1180  end subroutine prg_write_trajectoryandproperty
1181 
1182 
1190  subroutine prg_make_random_system(system,nats,seed,lx,ly,lz)
1191 
1192  implicit none
1193  integer :: i, nats, seed
1194  real(dp) :: lx, ly, lz, ran
1195  type(system_type), intent(out) :: system
1196 
1197  !Allocating the system
1198  system%nats = nats
1199  allocate(system%symbol(nats))
1200  allocate(system%atomic_number(nats))
1201  allocate(system%coordinate(3,nats))
1202  allocate(system%mass(nats))
1203 
1204  call random_seed(seed)
1205 
1206  do i=1,nats
1207 
1208  system%atomic_number(i) = 1
1209  call random_number(ran)
1210  system%coordinate(i,1) = ran*lx
1211  call random_number(ran)
1212  system%coordinate(i,2) = ran*ly
1213  call random_number(ran)
1214  system%coordinate(i,3) = ran*lz
1215 
1216  system%symbol(i) = element_symbol(system%atomic_number(i))
1217  system%mass(i) = element_mass(1)
1218 
1219  enddo
1220 
1221  allocate(system%lattice_vector(3,3))
1222  system%lattice_vector = 0.0_dp
1223  system%lattice_vector(1,1) = lx
1224  system%lattice_vector(2,2) = ly
1225  system%lattice_vector(3,3) = lz
1226 
1227  end subroutine prg_make_random_system
1228 
1236  subroutine prg_parameters_to_vectors(abc_angles,lattice_vector)
1237  implicit none
1238  real(dp) :: angle_alpha, angle_beta, angle_gamma, lattice_a
1239  real(dp) :: lattice_b, lattice_c, pi
1240  real(dp), intent(in) :: abc_angles(2,3)
1241  real(dp), intent(out) :: lattice_vector(3,3)
1242 
1243  pi = 3.14159265358979323846264338327950_dp
1244 
1245  lattice_a = abc_angles(1,1)
1246  lattice_b = abc_angles(1,2)
1247  lattice_c = abc_angles(1,3)
1248  angle_alpha = abc_angles(2,1)
1249  angle_beta = abc_angles(2,2)
1250  angle_gamma = abc_angles(2,3)
1251 
1252  angle_alpha = 2.0_dp*pi*angle_alpha/360.0_dp
1253  angle_beta = 2.0_dp*pi*angle_beta/360.0_dp
1254  angle_gamma = 2.0_dp*pi*angle_gamma/360.0_dp
1255 
1256  lattice_vector(1,1) = lattice_a
1257  lattice_vector(1,2)=0
1258  lattice_vector(1,3)=0
1259 
1260  lattice_vector(2,1)=lattice_b*cos(angle_gamma)
1261  lattice_vector(2,2)=lattice_b*sin(angle_gamma)
1262  lattice_vector(2,3)=0
1263 
1264  lattice_vector(3,1)=lattice_c*cos(angle_beta)
1265  lattice_vector(3,2)=lattice_c*( cos(angle_alpha)-cos(angle_gamma)* &
1266  cos(angle_beta) )/sin(angle_gamma)
1267  lattice_vector(3,3)=sqrt(lattice_c**2 - lattice_vector(3,1)**2 - lattice_vector(3,2)**2)
1268 
1269  end subroutine prg_parameters_to_vectors
1270 
1278  subroutine prg_vectors_to_parameters(lattice_vector,abc_angles)
1279  implicit none
1280  real(dp) :: angle_alpha, angle_beta, angle_gamma, lattice_a
1281  real(dp) :: lattice_b, lattice_c, pi
1282  real(dp), intent(in) :: lattice_vector(3,3)
1283  real(dp), intent(out) :: abc_angles(2,3)
1284 
1285  pi = 3.14159265358979323846264338327950_dp
1286 
1287  lattice_a = sqrt(lattice_vector(1,1)**2+lattice_vector(1,2)**2+lattice_vector(1,3)**2)
1288  lattice_b = sqrt(lattice_vector(2,1)**2+lattice_vector(2,2)**2+lattice_vector(2,3)**2)
1289  lattice_c = sqrt(lattice_vector(3,1)**2+lattice_vector(3,2)**2+lattice_vector(3,3)**2)
1290 
1291  angle_gamma = dot_product(lattice_vector(1,:),lattice_vector(2,:))/(lattice_a*lattice_b)
1292  angle_beta = dot_product(lattice_vector(1,:),lattice_vector(3,:))/(lattice_a*lattice_c)
1293  angle_alpha = dot_product(lattice_vector(2,:),lattice_vector(3,:))/(lattice_b*lattice_c)
1294 
1295 
1296  angle_alpha = 360.0_dp*acos(angle_alpha)/(2.0_dp*pi)
1297  angle_beta = 360.0_dp*acos(angle_beta)/(2.0_dp*pi)
1298  angle_gamma = 360.0_dp*acos(angle_gamma)/(2.0_dp*pi)
1299 
1300  abc_angles(1,1) = lattice_a
1301  abc_angles(1,2) = lattice_b
1302  abc_angles(1,3) = lattice_c
1303 
1304  abc_angles(2,1) = angle_alpha
1305  abc_angles(2,2) = angle_beta
1306  abc_angles(2,3) = angle_gamma
1307 
1308  end subroutine prg_vectors_to_parameters
1309 
1314  subroutine prg_get_origin(coords,origin)
1315  implicit none
1316  integer :: i
1317  real(dp) :: max_x, max_y, max_z, min_x
1318  real(dp) :: min_y, min_z
1319  real(dp), allocatable, intent(inout) :: origin(:)
1320  real(dp), intent(in) :: coords(:,:)
1321 
1322  if(.not.allocated(origin)) allocate(origin(3))
1323 
1324  max_x = -1.0d5 ; max_y = -1.0d5 ; max_z = -1.0d5 ;
1325  min_x = 1.0d5 ; min_y = 1.0d5 ; min_z = 1.0d5 ;
1326 
1327  ! Getting the system limits.
1328  do i=1,size(coords,dim=2)
1329  max_x = max(max_x,coords(1,i))
1330  min_x = min(min_x,coords(1,i))
1331  max_y = max(max_y,coords(2,i))
1332  min_y = min(min_y,coords(2,i))
1333  max_z = max(max_z,coords(3,i))
1334  min_z = min(min_z,coords(3,i))
1335  enddo
1336 
1337  origin(1) = min_x
1338  origin(2) = min_y
1339  origin(3) = min_z
1340 
1341  end subroutine prg_get_origin
1342 
1347  subroutine prg_get_distancematrix(coords,dmat)
1348  implicit none
1349  integer :: i,j,nats
1350  real(dp) :: norm2
1351  real(dp), intent(in) :: coords(:,:)
1352  real(dp), intent(out), allocatable :: dmat(:,:)
1353 
1354  nats = size(coords,dim=2)
1355  if(.not.allocated(dmat)) allocate(dmat(nats,nats))
1356 
1357  ! Getting the system limits.
1358  do i=1,nats
1359  do j=1,nats
1360  dmat(i,j) = prg_norm2(coords(:,i)-coords(:,j))
1361  enddo
1362  enddo
1363 
1364  end subroutine prg_get_distancematrix
1365 
1366 
1372  subroutine prg_translateandfoldtobox(coords,lattice_vectors,origin, verbose)
1373  implicit none
1374  integer :: i
1375  integer, intent(in), optional :: verbose
1376  real(dp) :: max_x, max_y, max_z, min_x
1377  real(dp) :: min_y, min_z
1378  real(dp), allocatable, intent(inout) :: origin(:),coords(:,:)
1379  real(dp), intent(in) :: lattice_vectors(:,:)
1380 
1381  if(present(verbose))then
1382  if(verbose >= 1)write(*,*)"In prg_translateandfoldtobox ..."
1383  endif
1384 
1385  max_x = -1.0d5 ; max_y = -1.0d5 ; max_z = -1.0d5 ;
1386  min_x = 1.0d5 ; min_y = 1.0d5 ; min_z = 1.0d5 ;
1387 
1388  ! Getting the system limits.
1389  do i=1,size(coords,dim=2)
1390  max_x = max(max_x,coords(1,i))
1391  min_x = min(min_x,coords(1,i))
1392  max_y = max(max_y,coords(2,i))
1393  min_y = min(min_y,coords(2,i))
1394  max_z = max(max_z,coords(3,i))
1395  min_z = min(min_z,coords(3,i))
1396  enddo
1397 
1398  if(.not.allocated(origin))then
1399  allocate(origin(3))
1400  origin(1) = min_x
1401  origin(2) = min_y
1402  origin(3) = min_z
1403  endif
1404 
1405  coords(1,:) = coords(1,:) - origin(1)
1406  coords(2,:) = coords(2,:) - origin(2)
1407  coords(3,:) = coords(3,:) - origin(3)
1408 
1409  coords(1,:) = mod(coords(1,:)*1000000,1000000*lattice_vectors(1,1))/1000000;
1410  coords(2,:) = mod(coords(2,:)*1000000,1000000*lattice_vectors(2,2))/1000000;
1411  coords(3,:) = mod(coords(3,:)*1000000,1000000*lattice_vectors(3,3))/1000000;
1412 
1413  origin = 0.0_dp
1414 
1415  !Get new origin.
1416  origin(1) = -1.0d-1 ; origin(2) = -1.0d-1; origin(3) = -1.0d-1
1417 
1418  end subroutine prg_translateandfoldtobox
1419 
1420 
1426  subroutine prg_centeratbox(coords,lattice_vectors,verbose)
1427  implicit none
1428  integer :: i, nats
1429  integer, intent(in), optional :: verbose
1430  real(dp) :: gc(3)
1431  real(dp), allocatable, intent(inout) :: coords(:,:)
1432  real(dp), intent(in) :: lattice_vectors(:,:)
1433 
1434  if(present(verbose))then
1435  if(verbose >= 1)write(*,*)"In prg_centeratbox ..."
1436  endif
1437 
1438  nats=size(coords,dim=2)
1439 
1440  gc= 0.0d0
1441 
1442  !$omp parallel do default(none) private(i) &
1443  !$omp shared(coords,nats) &
1444  !$omp reduction(+:gc)
1445  do i=1,nats
1446  gc=gc + coords(:,i)
1447  enddo
1448  !$omp end parallel do
1449 
1450  gc=gc/real(nats,dp)
1451 
1452  !$omp parallel do default(none) private(i) &
1453  !$omp shared(coords,lattice_vectors,nats, gc)
1454  do i=1,nats
1455  coords(1,i) = coords(1,i) + lattice_vectors(1,1)/2.0d0 - gc(1)
1456  coords(2,i) = coords(2,i) + lattice_vectors(2,2)/2.0d0 - gc(2)
1457  coords(3,i) = coords(3,i) + lattice_vectors(3,3)/2.0d0 - gc(3)
1458  enddo
1459  !$omp end parallel do
1460 
1461  end subroutine prg_centeratbox
1462 
1468  subroutine prg_wraparound(coords,lattice_vectors,index,verbose)
1469  implicit none
1470  integer :: i, nats
1471  integer, intent(in) :: index
1472  integer, intent(in), optional :: verbose
1473  real(dp), allocatable, intent(inout) :: coords(:,:)
1474  real(dp), allocatable :: origin(:)
1475  real(dp), intent(in) :: lattice_vectors(:,:)
1476 
1477  if(present(verbose))then
1478  if(verbose >= 1)write(*,*)"In prg_wraparound ..."
1479  endif
1480 
1481  if(.not.allocated(origin)) allocate(origin(3))
1482 
1483  nats=size(coords,dim=2)
1484 
1485  origin(1) = -coords(1,index) + lattice_vectors(1,1)/2.0_dp
1486  origin(2) = -coords(2,index) + lattice_vectors(2,2)/2.0_dp
1487  origin(3) = -coords(3,index) + lattice_vectors(3,3)/2.0_dp
1488 
1489  coords(1,:) = coords(1,:) + origin(1)
1490  coords(2,:) = coords(2,:) + origin(2)
1491  coords(3,:) = coords(3,:) + origin(3)
1492 
1493  !$omp parallel do default(none) private(i) &
1494  !$omp shared(coords,lattice_vectors,nats)
1495  do i=1,nats
1496  if(coords(1,i) > lattice_vectors(1,1))coords(1,i)=coords(1,i)-lattice_vectors(1,1)
1497  if(coords(2,i) > lattice_vectors(2,2))coords(2,i)=coords(2,i)-lattice_vectors(2,2)
1498  if(coords(3,i) > lattice_vectors(3,3))coords(3,i)=coords(3,i)-lattice_vectors(3,3)
1499  if(coords(1,i) < 0.0_dp)coords(1,i)=coords(1,i)+lattice_vectors(1,1)
1500  if(coords(2,i) < 0.0_dp)coords(2,i)=coords(2,i)+lattice_vectors(2,2)
1501  if(coords(3,i) < 0.0_dp)coords(3,i)=coords(3,i)+lattice_vectors(3,3)
1502  enddo
1503  !$end omp parallel do
1504 
1505  end subroutine prg_wraparound
1506 
1507 
1513  subroutine prg_translatetogeomcandfoldtobox(coords,lattice_vectors,origin)
1514  implicit none
1515  integer :: i
1516  real(dp), allocatable, intent(inout) :: origin(:),coords(:,:)
1517  real(dp), intent(in) :: lattice_vectors(:,:)
1518  real(dp) :: geomc(3)
1519 
1520  if(.not.allocated(origin)) allocate(origin(3))
1521 
1522  ! Getting the geometric center.
1523  do i=1,size(coords,dim=2)
1524  geomc = geomc + coords(:,i)
1525  enddo
1526 
1527  geomc = geomc/size(coords,dim=2)
1528 
1529  coords(1,:) = coords(1,:) - geomc(1)
1530  coords(2,:) = coords(2,:) - geomc(2)
1531  coords(3,:) = coords(3,:) - geomc(3)
1532 
1533  coords(1,:) = mod(coords(1,:)*1000000,1000000*lattice_vectors(1,1))/1000000;
1534  coords(2,:) = mod(coords(2,:)*1000000,1000000*lattice_vectors(2,2))/1000000;
1535  coords(3,:) = mod(coords(3,:)*1000000,1000000*lattice_vectors(3,3))/1000000;
1536 
1537  origin = 0.0_dp
1538 
1539  !Shift origin slightly
1540  origin(1) = -1.0d-1 ; origin(2) = -1.0d-1; origin(3) = -1.0d-1
1541 
1542  end subroutine prg_translatetogeomcandfoldtobox
1543 
1552  subroutine prg_replicate(coords,symbols,lattice_vectors,nx,ny,nz)
1553  implicit none
1554  integer :: i, j, k, l
1555  integer :: m, fnats, nats
1556  integer, intent(in) :: nx, ny, nz
1557  character(2), allocatable, intent(inout) :: symbols(:)
1558  character(2), allocatable :: rsymbols(:)
1559  real(dp), allocatable, intent(inout) :: coords(:,:)
1560  real(dp), intent(inout) :: lattice_vectors(:,:)
1561  real(dp), allocatable :: r(:,:)
1562 
1563  nats = size(coords,dim=2)
1564  fnats = nx*ny*nz*nats
1565 
1566  allocate(r(3,fnats))
1567  allocate(rsymbols(fnats))
1568 
1569  m = 0
1570  do i=1,nx
1571  do j=1,ny
1572  do k=1,nz
1573  do l=1,nats
1574  m=m+1
1575  r(1,m)=i*lattice_vectors(1,1) + j*lattice_vectors(1,2) + k*lattice_vectors(1,3)
1576  r(1,m)=r(1,m) + coords(1,l)
1577  r(2,m)=i*lattice_vectors(2,1) + j*lattice_vectors(2,2) + k*lattice_vectors(2,3)
1578  r(2,m)=r(2,m) + coords(2,l)
1579  r(3,m)=i*lattice_vectors(3,1) + j*lattice_vectors(3,2) + k*lattice_vectors(3,3)
1580  r(3,m)=r(3,m) + coords(3,l)
1581  rsymbols(m) = symbols(l)
1582  enddo
1583  enddo
1584  enddo
1585  enddo
1586 
1587  deallocate(coords,symbols)
1588  allocate(coords(3,fnats))
1589  allocate(symbols(fnats))
1590 
1591  lattice_vectors(1,:) = nx*lattice_vectors(1,:)
1592  lattice_vectors(2,:) = ny*lattice_vectors(2,:)
1593  lattice_vectors(3,:) = nz*lattice_vectors(3,:)
1594 
1595  coords = r
1596  symbols = rsymbols
1597  deallocate(r,rsymbols)
1598 
1599  end subroutine prg_replicate
1600 
1601 
1609  subroutine prg_replicate_system(sy,syf,nx,ny,nz)
1610  implicit none
1611  integer :: i, j, k, l
1612  integer :: ini, fin, cont, ii
1613  integer, intent(in) :: nx, ny, nz
1614  type(system_type), intent(inout) :: sy
1615  type(system_type) :: syf
1616 
1617  if(.not. allocated(sy%coordinate))then
1618  write(*,*)"ERROR: System does not have coordinate"
1619  stop
1620  endif
1621  if(.not. allocated(sy%symbol))then
1622  write(*,*)"ERROR: System does not have atomic symbols"
1623  stop
1624  endif
1625  syf%nats = sy%nats*nx*ny*nz
1626  allocate(syf%coordinate(3,syf%nats))
1627  allocate(syf%symbol(syf%nats))
1628  allocate(syf%resindex(syf%nats))
1629  allocate(syf%atomname(syf%nats))
1630  allocate(syf%atomic_number(syf%nats))
1631  allocate(syf%spindex(syf%nats))
1632 
1633  cont = 0
1634  do i = 1,nx
1635  do j = 1,ny
1636  do k = 1,nz
1637  do ii = 1,sy%nats
1638  cont = cont + 1
1639  syf%coordinate(1,cont) = sy%coordinate(1,ii) + &
1640  & sy%lattice_vector(1,1)*i + sy%lattice_vector(2,1)*i + sy%lattice_vector(3,1)*i
1641  syf%coordinate(2,cont) = sy%coordinate(2,ii) + &
1642  & sy%lattice_vector(1,2)*j + sy%lattice_vector(2,2)*j + sy%lattice_vector(3,2)*j
1643  syf%coordinate(3,cont) = sy%coordinate(3,ii) + &
1644  & sy%lattice_vector(1,3)*k + sy%lattice_vector(2,3)*k + sy%lattice_vector(3,3)*k
1645  syf%symbol(cont) = sy%symbol(ii)
1646  enddo
1647  enddo
1648  enddo
1649  enddo
1650 
1651  if(allocated(sy%resindex))then
1652  do i = 1,nx*ny*nz
1653  ini = (i - 1)*sy%nats + 1
1654  fin = i*sy%nats
1655  syf%resindex(ini:fin) = sy%resindex(:)
1656  enddo
1657  else
1658  syf%resindex = 1
1659  endif
1660 
1661  if(allocated(sy%atomname))then
1662  do i = 1,nx*ny*nz
1663  ini = (i - 1)*sy%nats + 1
1664  fin = i*sy%nats
1665  syf%atomname(ini:fin) = sy%atomname(:)
1666  enddo
1667  else
1668  syf%atomname = "At"
1669  endif
1670 
1671  if(allocated(sy%atomic_number))then
1672  do i = 1,nx*ny*nz
1673  ini = (i - 1)*sy%nats + 1
1674  fin = i*sy%nats
1675  syf%atomic_number(ini:fin) = sy%atomic_number(:)
1676  enddo
1677  endif
1678  if(allocated(sy%spindex))then
1679  do i = 1,nx*ny*nz
1680  ini = (i - 1)*sy%nats + 1
1681  fin = i*sy%nats
1682  syf%spindex(ini:fin) = sy%spindex(:)
1683  enddo
1684  endif
1685 
1686  allocate(syf%lattice_vector(3,3))
1687  syf%lattice_vector(1,:) = nx*sy%lattice_vector(1,:)
1688  syf%lattice_vector(2,:) = ny*sy%lattice_vector(2,:)
1689  syf%lattice_vector(3,:) = nz*sy%lattice_vector(3,:)
1690 
1691  if(allocated(sy%splist))then
1692  syf%nsp = sy%nsp
1693  allocate(syf%splist(syf%nsp))
1694  syf%splist = sy%splist
1695  endif
1696 
1697  end subroutine prg_replicate_system
1698 
1699 
1706  subroutine prg_cleanuprepeatedatoms(nats,coords,symbols,verbose)
1707  implicit none
1708  integer :: i, natstmp, l, count, j, count1
1709  integer, intent(inout) :: nats
1710  integer, intent(in), optional :: verbose
1711  real(dp) :: d
1712  real(dp), allocatable, intent(inout) :: coords(:,:)
1713  real(dp), allocatable :: coordstmp(:,:)
1714  character(len=*), allocatable, intent(inout) :: symbols(:)
1715  character(2), allocatable :: symbolstmp(:)
1716 
1717  if(present(verbose))then
1718  if(verbose >= 1)write(*,*)"In prg_cleanuprepeatedatoms ..."
1719  endif
1720 
1721  natstmp=nats
1722 
1723  allocate(coordstmp(3,nats))
1724  allocate(symbolstmp(nats))
1725 
1726  coordstmp(1,1)=coords(1,1)
1727  coordstmp(2,1)=coords(2,1)
1728  coordstmp(3,1)=coords(3,1)
1729 
1730  l=1
1731  symbolstmp(1)=symbols(1)
1732 
1733  do j=1,nats
1734  count=0
1735  count1=0
1736  do i=1,l
1737  d=sqrt((coords(1,j)-coordstmp(1,i))**2+(coords(2,j)-coordstmp(2,i))**2+&
1738  &(coords(3,j)-coordstmp(3,i))**2)
1739  if(d.gt.0.1d0)then
1740  count=count + 1
1741  else
1742  count1=count1 + 1
1743  endif
1744  enddo
1745 
1746  if(count.eq.l)then ! New atom in collection
1747  if(count1 >= 1)then
1748  write(*,*)"Problem at iteration",j
1749  stop
1750  endif
1751 
1752  l = l+1
1753  coordstmp(1,l)=coords(1,j)
1754  coordstmp(2,l)=coords(2,j)
1755  coordstmp(3,l)=coords(3,j)
1756  symbolstmp(l)=symbols(j)
1757  endif
1758  enddo
1759 
1760  deallocate(symbols)
1761  deallocate(coords)
1762  nats = l
1763  allocate(symbols(nats))
1764  allocate(coords(3,nats))
1765 
1766  do i=1,nats
1767  symbols(i) = symbolstmp(i)
1768  coords(:,i) = coordstmp(:,i)
1769  enddo
1770 
1771  deallocate(coordstmp)
1772  deallocate(symbolstmp)
1773 
1774  end subroutine prg_cleanuprepeatedatoms
1775 
1787  subroutine prg_get_recip_vects(lattice_vectors,recip_vectors,volr,volk)
1788  implicit none
1789  real(dp) :: a1xa2(3), a2xa3(3), a3xa1(3), b2xb3(3)
1790  real(dp) :: pi
1791  real(dp), allocatable, intent(inout) :: recip_vectors(:,:)
1792  real(dp), intent(in) :: lattice_vectors(:,:)
1793  real(dp), intent(inout) :: volk, volr
1794 
1795  if(.not.allocated(recip_vectors))allocate(recip_vectors(3,3))
1796  volr=0.0_dp; volk=0.0_dp
1797 
1798  pi = 3.14159265358979323846264338327950_dp
1799 
1800  a1xa2(1) = lattice_vectors(1,2)*lattice_vectors(2,3) - lattice_vectors(1,3)*lattice_vectors(2,2)
1801  a1xa2(2) = -lattice_vectors(1,1)*lattice_vectors(2,3) + lattice_vectors(1,3)*lattice_vectors(2,1)
1802  a1xa2(3) = lattice_vectors(1,1)*lattice_vectors(2,2) - lattice_vectors(1,2)*lattice_vectors(2,1)
1803 
1804  a2xa3(1) = lattice_vectors(2,2)*lattice_vectors(3,3) - lattice_vectors(2,3)*lattice_vectors(3,2)
1805  a2xa3(2) = -lattice_vectors(2,1)*lattice_vectors(3,3) + lattice_vectors(2,3)*lattice_vectors(3,1)
1806  a2xa3(3) = lattice_vectors(2,1)*lattice_vectors(3,2) - lattice_vectors(2,2)*lattice_vectors(3,1)
1807 
1808  a3xa1(1) = lattice_vectors(3,2)*lattice_vectors(1,3) - lattice_vectors(3,3)*lattice_vectors(1,2)
1809  a3xa1(2) = -lattice_vectors(3,1)*lattice_vectors(1,3) + lattice_vectors(3,3)*lattice_vectors(1,1)
1810  a3xa1(3) = lattice_vectors(3,1)*lattice_vectors(1,2) - lattice_vectors(3,2)*lattice_vectors(1,1)
1811 
1812  !Get the volume of the cell
1813  volr = lattice_vectors(1,1)*a2xa3(1)+ lattice_vectors(1,2)*a2xa3(2)+lattice_vectors(1,3)*a2xa3(3)
1814 
1815  !Get the reciprocal vectors
1816  recip_vectors(1,:) = 2.0_dp*pi*a2xa3(:)/volr
1817  recip_vectors(2,:) = 2.0_dp*pi*a3xa1(:)/volr
1818  recip_vectors(3,:) = 2.0_dp*pi*a1xa2(:)/volr
1819 
1820  b2xb3(1) = recip_vectors(2,2)*recip_vectors(3,3) - recip_vectors(2,3)*recip_vectors(3,2)
1821  b2xb3(2) = -recip_vectors(2,1)*recip_vectors(3,3) + recip_vectors(2,3)*recip_vectors(3,1)
1822  b2xb3(3) = recip_vectors(2,1)*recip_vectors(3,2) - recip_vectors(2,2)*recip_vectors(3,1)
1823 
1824  volk = recip_vectors(1,1)*b2xb3(1)+ recip_vectors(1,2)*b2xb3(2)+recip_vectors(1,3)*b2xb3(3)
1825 
1826  end subroutine prg_get_recip_vects
1827 
1836  subroutine prg_get_dihedral(coords,id1,id2,id3,id4,dihedral)
1837 
1838  real(dp) :: mv1, mv2, v1(3), v2(3)
1839  real(dp) :: dotprod, cosdir, v2xv20(3), v1xv10(3)
1840  real(dp) :: v10(3),v20(3), cprod(3), normcprod, sindir
1841  real(dp), intent(in) :: coords(:,:)
1842  real(dp), intent(out) :: dihedral
1843  integer :: i
1844  integer, intent(in) :: id1,id2,id3,id4
1845  character(2) :: index1, index2, index3, index4
1846 
1847  v1=coords(:,id4) - coords(:,id3)
1848  v10=coords(:,id2) - coords(:,id3)
1849  v2=coords(:,id1) - coords(:,id2)
1850  v20=coords(:,id3) - coords(:,id2)
1851 
1852  v1xv10(1)=v1(2)*v10(3)-v1(3)*v10(2)
1853  v1xv10(2)=-(v1(1)*v10(3)-v1(3)*v10(1))
1854  v1xv10(3)=v1(1)*v10(2)-v1(2)*v10(1)
1855 
1856  v2xv20(1)=v2(2)*v20(3)-v2(3)*v20(2)
1857  v2xv20(2)=-(v2(1)*v20(3)-v2(3)*v20(1))
1858  v2xv20(3)=v2(1)*v20(2)-v2(2)*v20(1)
1859 
1860  dotprod = v1xv10(1)*v2xv20(1) + v1xv10(2)*v2xv20(2) + v1xv10(3)*v2xv20(3)
1861 
1862  cprod(1)=v1xv10(2)*v2xv20(3)-v1xv10(3)*v2xv20(2)
1863  cprod(2)=-(v1xv10(1)*v2xv20(3)-v1xv10(3)*v2xv20(1))
1864  cprod(3)=v1xv10(1)*v2xv20(2)-v1xv10(2)*v2xv20(1)
1865 
1866  normcprod=sqrt(cprod(1)*cprod(1) + cprod(2)*cprod(2) + cprod(3)*cprod(3))
1867 
1868  mv1= sqrt(v1xv10(1)*v1xv10(1) + v1xv10(2)*v1xv10(2) + v1xv10(3)*v1xv10(3))
1869  mv2= sqrt(v2xv20(1)*v2xv20(1) + v2xv20(2)*v2xv20(2) + v2xv20(3)*v2xv20(3))
1870 
1871  cosdir = dotprod/(mv1*mv2)
1872  sindir = normcprod/(mv1*mv2)
1873 
1874  sindir = cprod(3)
1875  dihedral=sign(1.0d0,sindir)*acos(-cosdir)
1876  dihedral=-360*dihedral/(2.0*3.14159265359)
1877  if(dihedral < 0)dihedral = 360 + dihedral
1878 
1879  end subroutine prg_get_dihedral
1880 
1891  subroutine prg_get_covgraph(sy,nnStruct,nrnnstruct,bml_type,factor,gcov_bml,mdimin,verbose)
1892  implicit none
1893  character(20), intent(in) :: bml_type
1894  integer :: i, j, jj, mymdim
1895  integer, intent(in) :: mdimin
1896  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:)
1897  integer, optional, intent(in) :: verbose
1898  real(dp) :: d, dvdw, factor, ra(3),rb(3),rab(3)
1899  real(dp) :: lx, ly, lz
1900  type(bml_matrix_t), intent(inout) :: gcov_bml
1901  type(system_type), intent(in) :: sy
1902  logical(1), allocatable :: ispresent(:)
1903  real(kind(1.0)), allocatable :: row(:)
1904 
1905  if(bml_get_n(gcov_bml).gt.0) call bml_deallocate(gcov_bml)
1906 
1907  if(mdimin > 0)then
1908  mymdim = mdimin
1909  else
1910  mymdim = sy%nats
1911  endif
1912 
1913  call bml_zero_matrix(bml_type,bml_element_real,kind(1.0),sy%nats,mymdim,gcov_bml)
1914  !call bml_zero_matrix(bml_type,bml_element_real,dp,sy%nats,mymdim,gcov_bml)
1915 
1916  lx = sy%lattice_vector(1,1)
1917  ly = sy%lattice_vector(2,2)
1918  lz = sy%lattice_vector(3,3)
1919 
1920  if(present(verbose))then
1921  if(verbose >= 1)then
1922  write(*,*)" "
1923  write(*,*)"Building covalency graph ..."
1924  write(*,*)" "
1925  endif
1926  endif
1927 
1928  allocate(ispresent(sy%nats))
1929 
1930  do i=1,sy%nats
1931  ra = sy%coordinate(:,i)
1932  ispresent = .false.
1933  do j=1,nrnnstruct(i)
1934 
1935  if(nrnnstruct(i) <= 1)write(*,*)"WARNING: Atom" ,i,"is desconnected"
1936 
1937  jj = nnstruct(j,i)
1938  if(.not.ispresent(jj))then
1939  rb = sy%coordinate(:,jj)
1940 
1941  rab(1) = modulo((rb(1) - ra(1) + lx/2.0_dp),lx) - lx/2.0_dp
1942  rab(2) = modulo((rb(2) - ra(2) + ly/2.0_dp),ly) - ly/2.0_dp
1943  rab(3) = modulo((rb(3) - ra(3) + lz/2.0_dp),lz) - lz/2.0_dp
1944 
1945  d = prg_norm2(rab)
1946 
1947  if(d == 0.0_dp .and. i .ne. jj)write(*,*)"WARNING: Atom" ,i,"and atom",jj,&
1948  "are on top of each other"
1949  dvdw = factor*(element_covr(sy%atomic_number(i)) + element_covr(sy%atomic_number(jj)))
1950  if(d < dvdw .and. d > 0.0_dp)then
1951  ispresent(jj) = .true.
1952  ! call bml_set_element_new(gcov_bml,i,jj,1.0_dp)
1953  ! call bml_set_element_new(gcov_bml,jj,i,1.0_dp)
1954  call bml_set_element_new(gcov_bml,i,jj,1.0)
1955  ! call bml_set_element_new(gcov_bml,jj,i,1.0)
1956  endif
1957  endif
1958  enddo
1959 
1960  ! call bml_set_element_new(gcov_bml,i,i,1.0_dp)
1961  enddo
1962 
1963  deallocate(ispresent)
1964 
1965  end subroutine prg_get_covgraph
1966 
1967 
1968  subroutine prg_get_covgraph_int(sy,nnStructMindist,nnStruct,nrnnstruct,bml_type,factor,gcov_bml,mdimin,verbose)
1969  implicit none
1970  character(20), intent(in) :: bml_type
1971  integer :: i, j, jj, mymdim
1972  integer, intent(in) :: mdimin
1973  integer, intent(in) :: nnStruct(:,:), nrnnstruct(:)
1974  integer, optional, intent(in) :: verbose
1975  real(dp) :: d, dvdw, factor
1976  real(dp), intent(in) :: nnStructMindist(:,:)
1977  type(bml_matrix_t), intent(inout) :: gcov_bml
1978  type(system_type), intent(in) :: sy
1979 
1980  if(bml_get_n(gcov_bml).gt.0) call bml_deallocate(gcov_bml)
1981 
1982  if(mdimin > 0)then
1983  mymdim = mdimin
1984  else
1985  mymdim = sy%nats
1986  endif
1987 
1988  call bml_zero_matrix(bml_type,bml_element_real,dp,sy%nats,mymdim,gcov_bml)
1989 
1990  if(present(verbose))then
1991  if(verbose >= 1)then
1992  write(*,*)" "
1993  write(*,*)"Building covalency graph ..."
1994  write(*,*)" "
1995  endif
1996  endif
1997 
1998  do i=1,sy%nats
1999  do j=1,nrnnstruct(i)
2000  if(nrnnstruct(i) <= 1)write(*,*)"WARNING: Atom" ,i,"is desconnected"
2001  jj = nnstruct(j,i)
2002  d = nnstructmindist(j,i)
2003  if(d == 0.0_dp .and. i .ne. jj)write(*,*)"WARNING: Atom" ,i,"and atom",j,&
2004  "are on top of each other"
2005  dvdw = factor*(element_covr(sy%atomic_number(i)) + element_covr(sy%atomic_number(jj)))
2006  if(d < dvdw .and. d > 0.0_dp)then
2007  call bml_set_element_new(gcov_bml,i,jj,1.0_dp)
2008  call bml_set_element_new(gcov_bml,jj,i,1.0_dp)
2009  endif
2010  enddo
2011  call bml_set_element_new(gcov_bml,i,i,1.0_dp)
2012  enddo
2013 
2014  end subroutine prg_get_covgraph_int
2015 
2016 
2027  subroutine prg_get_covgraph_h(sy,nnStruct,nrnnstruct,rcut,&
2028  graph_h,mdimin,verbose)
2029  implicit none
2030  integer :: i, j, jj, ncount, mymdim
2031  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:), mdimin
2032  integer, optional, intent(in) :: verbose
2033  real(dp) :: d, dvdw, ra(3),rb(3),rab(3)
2034  real(dp) :: lx, ly, lz
2035  real(dp), intent(in) :: rcut
2036  integer, allocatable, intent(inout) :: graph_h(:,:)
2037  type(system_type), intent(in) :: sy
2038 
2039 
2040  if(mdimin > 0)then
2041  mymdim = mdimin
2042  else
2043  mymdim = sy%nats
2044  endif
2045 
2046  if(.not.allocated(graph_h)) allocate(graph_h(mymdim,sy%nats))
2047 
2048  if(present(verbose))then
2049  if(verbose >= 1)then
2050  write(*,*)" "
2051  write(*,*)"Building H covalency graph ..."
2052  write(*,*)" "
2053  endif
2054  endif
2055 
2056  graph_h = 0
2057 
2058  lx = sy%lattice_vector(1,1)
2059  ly = sy%lattice_vector(2,2)
2060  lz = sy%lattice_vector(3,3)
2061 
2062  !$omp parallel do default(none) private(i) &
2063  !$omp private(ncount,j,jj,ra,rb,rab,d) &
2064  !$omp shared(sy,nrnnstruct,nnStruct,Lx,Ly,Lz,graph_h,rcut)
2065  do i=1,sy%nats
2066  ncount = 0
2067  do j=1,nrnnstruct(i)
2068  jj = nnstruct(j,i)
2069  ra=sy%coordinate(:,i)
2070  rb=sy%coordinate(:,jj)
2071  rab(1) = modulo((rb(1) - ra(1) + lx/2.0_dp),lx) - lx/2.0_dp
2072  rab(2) = modulo((rb(2) - ra(2) + ly/2.0_dp),ly) - ly/2.0_dp
2073  rab(3) = modulo((rb(3) - ra(3) + lz/2.0_dp),lz) - lz/2.0_dp
2074 
2075  d = prg_norm2(rab)
2076 
2077  if(d.lt.rcut.and.d.gt.0.0_dp)then
2078  ncount=ncount+1
2079  graph_h(ncount,i) = jj
2080  endif
2081  enddo
2082  enddo
2083  !$omp end parallel do
2084 
2085  end subroutine prg_get_covgraph_h
2086 
2097  subroutine prg_get_subsystem(sy,lsize,indices,sbsy,verbose)
2098  implicit none
2099  integer :: i, nsptmp, prev
2100  integer, intent(in) :: indices(:), lsize
2101  integer, optional, intent(in) :: verbose
2102  type(system_type), intent(in) :: sy
2103  type(system_type), intent(inout) :: sbsy
2104 
2105  if(present(verbose))then
2106  if(verbose >= 1)then
2107  write(*,*)" "
2108  write(*,*)"Extracting subsystem ..."
2109  write(*,*)" "
2110  endif
2111  endif
2112 
2113  sbsy%nats = lsize
2114 
2115  if(allocated(sbsy%symbol))then
2116  deallocate(sbsy%symbol)
2117  deallocate(sbsy%coordinate)
2118  deallocate(sbsy%atomic_number)
2119  deallocate(sbsy%velocity)
2120  deallocate(sbsy%force)
2121  deallocate(sbsy%net_charge)
2122  deallocate(sbsy%mass)
2123  deallocate(sbsy%lattice_vector)
2124  deallocate(sbsy%spindex)
2125  endif
2126 
2127  allocate(sbsy%symbol(sbsy%nats))
2128  allocate(sbsy%coordinate(3,sbsy%nats))
2129  allocate(sbsy%atomic_number(sbsy%nats))
2130  allocate(sbsy%velocity(3,sbsy%nats))
2131  allocate(sbsy%force(3,sbsy%nats))
2132  allocate(sbsy%net_charge(sbsy%nats))
2133  allocate(sbsy%mass(sbsy%nats))
2134  allocate(sbsy%lattice_vector(3,3))
2135  allocate(sbsy%spindex(sbsy%nats))
2136 
2137  ! allocate(sbsy%recip_vector(3,3))
2138 
2139  sbsy%lattice_vector = sy%lattice_vector
2140  sbsy%volr = sy%volr
2141  sbsy%volk = sy%volk
2142  nsptmp = 0
2143  do i=1,sbsy%nats
2144  sbsy%symbol(i) = sy%symbol(indices(i)+1) !Indices from the graph partition start from 0
2145  sbsy%coordinate(:,i) = sy%coordinate(:,indices(i)+1)
2146  sbsy%atomic_number(i) = sy%atomic_number(indices(i)+1)
2147  sbsy%mass(i) = sy%mass(indices(i)+1)
2148  sbsy%spindex(i) = sy%spindex(indices(i)+1)
2149  if(sbsy%spindex(i).gt.nsptmp)then
2150  nsptmp = nsptmp + 1
2151  endif
2152  enddo
2153 
2154  if(nsptmp.lt.sy%nsp)then
2155  write(*,*)"WARNING: nsp = ",nsptmp,sy%nsp
2156  write(*,*)"The subsystem contains less species that the system ..."
2157  write(*,*)"A generalization to parts where subsystem contains"
2158  write(*,*)"less species that the system needs to be added ..."
2159  write(*,*)"See prg_get_subsystem in prg_system_mod"
2160  ! stop
2161  endif
2162 
2163  sbsy%nsp = sy%nsp
2164 
2165  if(allocated(sbsy%spatnum))then
2166  deallocate(sbsy%spatnum)
2167  deallocate(sbsy%spmass)
2168  endif
2169 
2170  allocate(sbsy%spatnum(sbsy%nsp))
2171  allocate(sbsy%spmass(sbsy%nsp))
2172 
2173  sbsy%spatnum = sy%spatnum
2174  sbsy%spmass = sy%spmass
2175 
2176  end subroutine prg_get_subsystem
2177 
2178 
2185  subroutine prg_destroy_subsystems(sbsy,verbose)
2186  implicit none
2187  type(system_type), intent(inout) :: sbsy
2188  integer, optional, intent(in) :: verbose
2189  integer :: nparts
2190 
2191  if(present(verbose))then
2192  if(verbose >= 1)then
2193  write(*,*)" "
2194  write(*,*)"Deallocating the multiple subsystem ..."
2195  write(*,*)" "
2196  endif
2197  endif
2198 
2199 
2200  if(allocated(sbsy%symbol))then
2201  deallocate(sbsy%symbol)
2202  deallocate(sbsy%coordinate)
2203  deallocate(sbsy%atomic_number)
2204  deallocate(sbsy%velocity)
2205  deallocate(sbsy%force)
2206  deallocate(sbsy%net_charge)
2207  deallocate(sbsy%mass)
2208  deallocate(sbsy%spindex)
2209  endif
2210 
2211  if(allocated(sbsy%lattice_vector))deallocate(sbsy%lattice_vector)
2212  if(allocated(sbsy%resindex))deallocate(sbsy%resindex)
2213 
2214  if(allocated(sbsy%spatnum))then
2215  deallocate(sbsy%spatnum)
2216  deallocate(sbsy%spmass)
2217  endif
2218 
2219  if(bml_get_n(sbsy%estr%ham) > 0) call bml_deallocate(sbsy%estr%ham)
2220  if(bml_get_n(sbsy%estr%ham0) > 0) call bml_deallocate(sbsy%estr%ham0)
2221  if(bml_get_n(sbsy%estr%oham) > 0) call bml_deallocate(sbsy%estr%oham)
2222  if(bml_get_n(sbsy%estr%over) > 0) call bml_deallocate(sbsy%estr%over)
2223  if(bml_get_n(sbsy%estr%rho) > 0) call bml_deallocate(sbsy%estr%rho)
2224  if(bml_get_n(sbsy%estr%orho) > 0) call bml_deallocate(sbsy%estr%orho)
2225  if(bml_get_n(sbsy%estr%zmat) > 0) call bml_deallocate(sbsy%estr%zmat)
2226 
2227 
2228  if(allocated(sbsy%estr%coul_pot_r)) deallocate(sbsy%estr%coul_pot_r)
2229  if(allocated(sbsy%estr%coul_pot_k)) deallocate(sbsy%estr%coul_pot_k)
2230  if(allocated(sbsy%estr%skforce)) deallocate(sbsy%estr%skforce)
2231  if(allocated(sbsy%estr%fpul)) deallocate(sbsy%estr%fpul)
2232  if(allocated(sbsy%estr%fscoul)) deallocate(sbsy%estr%fscoul)
2233  if(allocated(sbsy%estr%hindex)) deallocate(sbsy%estr%hindex)
2234 
2235 
2236  end subroutine prg_destroy_subsystems
2237 
2249  subroutine prg_molpartition(sy,npart,nnStructMindist,nnStruct,nrnnstruct,hetatm,gp,verbose)
2250  implicit none
2251  character(2), intent(in) :: hetatm
2252  integer :: cnt, i, j, jj
2253  integer :: maxparts, nmax
2254  integer, allocatable :: core_count(:), cores(:,:), partof(:)
2255  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:)
2256  integer, intent(inout) :: npart
2257  integer, optional, intent(inout) :: verbose
2258  logical, allocatable :: inpart(:)
2259  real(dp) :: d, dvdw
2260  real(dp), intent(in) :: nnstructmindist(:,:)
2261  type(graph_partitioning_t), intent(inout) :: gp
2262  type(system_type), intent(in) :: sy
2263 
2264  if(present(verbose))then
2265  if(verbose >= 1)then
2266  write(*,*)" "; write(*,*)"Partitioning by molecule ..."; write(*,*)" "
2267  endif
2268  endif
2269 
2270  if(allocated(gp%nnodesInPart))then
2271  call prg_destroygraphpartitioning(gp)
2272  endif
2273 
2274  nmax = sy%nats
2275  maxparts = sy%nats
2276 
2277  allocate(core_count(maxparts))
2278 
2279  allocate(cores(nmax,maxparts))
2280  allocate(inpart(sy%nats))
2281  allocate(partof(sy%nats))
2282 
2283  cores = -1
2284  inpart = .false.
2285  partof=-1
2286  core_count = 0
2287  npart = 0
2288  do i=1,sy%nats
2289  if(.not.inpart(i).and.trim(sy%symbol(i)) == trim(hetatm))then
2290  npart = npart + 1
2291  core_count(npart) = core_count(npart) + 1
2292  cores(core_count(npart),npart) = i
2293  partof(i) = npart
2294  do j=1,nrnnstruct(i)
2295  jj = nnstruct(j,i)
2296  d = nnstructmindist(j,i)
2297  dvdw = 1.3_dp
2298  if(d.lt.dvdw.and.d.gt.0.001_dp)then
2299  if(.not.inpart(jj))then
2300  core_count(npart) = core_count(npart) + 1
2301  cores(core_count(npart),npart) = jj
2302  inpart(jj) = .true.
2303  partof(jj) = npart
2304  endif
2305  endif
2306  enddo
2307  endif
2308  enddo
2309 
2310  call prg_initgraphpartitioning(gp, "molecules", npart, sy%nats, sy% nats)
2311 
2312  ! Initialize and fill up subgraph structure
2313  ! Assign node ids (mapped to orbitals as rows) to each node in each
2314  do i = 1, npart
2315  gp%nnodesInPartAll(i) = core_count(i)
2316  call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
2317  allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
2318  gp%nnodesInPart(i) = core_count(i)
2319  enddo
2320 
2321  !Assign node ids to sgraph
2322  do i=1, npart
2323  do j=1,core_count(i)
2324  gp%sgraph(i)%nodeInPart(j) = cores(j,i) - 1
2325  end do
2326  enddo
2327 
2328  end subroutine prg_molpartition
2329 
2338  subroutine prg_get_partial_atomgraph(rho_bml,hindex,gch_bml,threshold,verbose)
2339  implicit none
2340  character(20) :: bml_type
2341  integer :: i, ii, j, jj
2342  integer :: nch, norbs
2343  integer, intent(in) :: hindex(:,:)
2344  integer, optional, intent(in) :: verbose
2345  logical :: connection
2346  logical, allocatable :: iconnectedtoj(:)
2347  real(dp), allocatable :: row(:), rowat(:)
2348  real(dp), intent(in) :: threshold
2349  type(bml_matrix_t), intent(in) :: rho_bml
2350  type(bml_matrix_t), intent(inout) :: gch_bml
2351 
2352  norbs = bml_get_n(rho_bml)
2353  bml_type = bml_get_type(rho_bml)
2354  nch = size(hindex,dim=2)
2355  allocate(rowat(nch))
2356  allocate(row(norbs))
2357  allocate(iconnectedtoj(nch))
2358 
2359  if(bml_get_n(gch_bml) > 0) call bml_deallocate(gch_bml)
2360  call bml_zero_matrix(bml_type,bml_element_real,dp,nch,nch,gch_bml)
2361 
2362  do i = 1, nch
2363  rowat = 0.0_dp
2364  do ii = hindex(1,i),hindex(2,i) ! i,j block
2365  call bml_get_row(rho_bml,ii,row)
2366  iconnectedtoj = .false.
2367  do j = 1, nch
2368  connection = .false.
2369 
2370  if(.not.iconnectedtoj(j).and.rowat(j) == 0)then
2371  do jj = hindex(1,j),hindex(2,j)
2372  if(abs(row(jj)) > threshold)then
2373  connection = .true. !We exit if there is a connection.
2374  exit
2375  endif
2376  enddo
2377  if(connection) iconnectedtoj(j) = .true.
2378  endif
2379  if(connection) rowat(j) = 1.0_dp
2380  enddo
2381 
2382  enddo
2383  call bml_set_row(gch_bml,i,rowat)
2384  enddo
2385 
2386  deallocate(rowat)
2387  deallocate(row)
2388  deallocate(iconnectedtoj)
2389 
2390  end subroutine prg_get_partial_atomgraph
2391 
2392 
2404  subroutine prg_collect_graph_p(rho_bml,nc,nats,hindex,chindex,graph_p,threshold,mdimin,verbose)
2405  implicit none
2406  character(20) :: bml_type
2407  integer :: i, ifull, ii, j
2408  integer :: jfull, jj, ncounti
2409  integer :: norbs, mymdim
2410  logical(1), allocatable :: rowatfull(:)
2411  integer, allocatable, intent(inout) :: graph_p(:,:)
2412  integer, intent(in) :: chindex(:), hindex(:,:), nats, nc
2413  integer :: nch
2414  integer, intent(in) :: mdimin
2415  integer, optional, intent(in) :: verbose
2416  logical :: connection
2417  logical, allocatable :: iconnectedtoj(:)
2418  real(dp), allocatable :: row(:)
2419  real(dp), intent(in) :: threshold
2420  type(bml_matrix_t), intent(in) :: rho_bml
2421 
2422  norbs = bml_get_n(rho_bml)
2423  bml_type = bml_get_type(rho_bml)
2424  nch = size(hindex,dim=2)
2425  allocate(rowatfull(nats))
2426  allocate(row(norbs))
2427  allocate(iconnectedtoj(nch))
2428 
2429  if(mdimin > 0)then
2430  mymdim = mdimin
2431  else
2432  mymdim = nats
2433  endif
2434 
2435  if(.not.allocated(graph_p))then
2436  allocate(graph_p(mymdim,nats))
2437  write(*,*)"PRG_COLLECT_GRAPH_P: Allocated graph_p"
2438  graph_p = 0
2439  endif
2440 
2441  !$omp parallel do default(none) private(i) &
2442  !$omp private(ncounti,rowatfull,ii,j,jfull,ifull,iconnectedtoj) &
2443  !$omp private(row,connection) &
2444  !$omp shared(graph_p,nc,chindex,hindex,rho_bml,nch,threshold,nats)
2445  do i = 1, nc
2446 
2447  ifull = chindex(i) + 1 !Map it to the full system
2448  if(ifull.gt.nats)write(*,*)"PRG_COLLECT_GRAPH_P: ifull > nats"
2449  rowatfull = .false.
2450  ii=1
2451  ncounti = 0
2452  do while (graph_p(ii,ifull).ne.0) !Unfolding the connections of ifull
2453  rowatfull(graph_p(ii,ifull)) = .true.
2454  ncounti = ncounti + 1
2455  ii = ii+1
2456  enddo
2457  do ii = hindex(1,i),hindex(2,i) ! i,j block
2458  call bml_get_row(rho_bml,ii,row)
2459  iconnectedtoj = .false.
2460  !
2461  do j = 1, nch
2462  jfull = chindex(j) + 1 !Cause the count starts form 0
2463  if(jfull.gt.nats)write(*,*)"PRG_COLLECT_GRAPH_P: jfull > nats"
2464  connection = .false.
2465  if(.not.iconnectedtoj(j).and..not.rowatfull(jfull))then
2466  do jj = hindex(1,j),hindex(2,j)
2467  if(abs(row(jj)) > threshold)then
2468  connection = .true. !We exit if there is a connection.
2469  exit
2470  endif
2471  enddo
2472  if(connection) iconnectedtoj(j) = .true.
2473  endif
2474  if(connection)then !I conected to j >>> ifull connected to jfull
2475  !Map back to the full system graph.
2476  rowatfull(jfull) = .true.
2477  ncounti = ncounti + 1
2478  graph_p(ncounti,ifull) = jfull
2479  endif
2480  enddo
2481  enddo
2482  enddo
2483  !$omp end parallel do
2484 
2485  deallocate(rowatfull)
2486  deallocate(iconnectedtoj)
2487  deallocate(row)
2488 
2489  end subroutine prg_collect_graph_p
2490 
2502  subroutine prg_collect_extended_graph_p(rho_bml,nc,nats,hindex,chindex,graph_p,threshold,mdimin,alpha,coords,coordsall,latticevectors,verbose)
2503  implicit none
2504  character(20) :: bml_type
2505  integer :: i, ifull, ii, j, k
2506  integer :: jfull, jj, nch, ncounti
2507  integer :: norbs, mymdim,chindex_size
2508  logical(1), allocatable :: rowatfull(:)
2509  integer, allocatable, intent(inout) :: graph_p(:,:)
2510  integer, intent(in) :: chindex(:), hindex(:,:), nats, nc
2511  integer, intent(in) :: mdimin
2512  integer, optional, intent(in) :: verbose
2513  logical :: connection
2514  logical, allocatable :: iconnectedtoj(:)
2515 #ifdef USE_SINGLE
2516  real(4), allocatable :: row(:),extmat(:,:),dvec(:,:),dr2(:)
2517  real(4), allocatable :: rho(:,:), rho_red(:,:), rhoext(:,:)
2518 #else
2519  real(dp), allocatable :: row(:),extmat(:,:),dvec(:,:),dr2(:)
2520  real(dp), allocatable :: rho(:,:), rho_red(:,:), rhoext(:,:)
2521 #endif
2522  integer, allocatable :: graph_core(:,:)
2523  real(dp), intent(in) :: alpha,threshold
2524  real(dp), allocatable, intent(in) :: latticevectors(:,:)
2525  real(dp), allocatable, intent(in) :: coords(:,:),coordsall(:,:)
2526  real(dp) :: lx, ly, lz, dvx, dvy, dvz, val
2527  type(bml_matrix_t), intent(inout) :: rho_bml
2528 #ifdef USE_OFFLOAD
2529  type(c_ptr) :: rho_bml_c_ptr
2530  integer :: ld
2531  real(c_double), pointer :: rho_bml_ptr(:,:)
2532  logical(1), allocatable :: graphed(:,:)
2533 #endif
2534 
2535  chindex_size = size(chindex)
2536  norbs = bml_get_n(rho_bml)
2537  bml_type = bml_get_type(rho_bml)
2538  nch = size(hindex,dim=2)
2539  allocate(rowatfull(nats))
2540  allocate(row(norbs))
2541  allocate(iconnectedtoj(nats))
2542  allocate(dvec(3,nats))
2543  allocate(dr2(nats))
2544  allocate(extmat(nats,nch))
2545  allocate(rho_red(nch,nc))
2546 
2547  if(mdimin > 0)then
2548  if(mdimin>nats)then
2549  mymdim = nats
2550  else
2551  mymdim = mdimin
2552  endif
2553  else
2554  mymdim = nats
2555  endif
2556 
2557 
2558 
2559  if(.not.allocated(graph_p))then
2560  write(*,*)"PRG_COLLECT_GRAPH_P: graph_p not allocated on entry, so allocating"
2561  allocate(graph_p(mymdim,nats))
2562  graph_p = 0
2563  else
2564  if((size(graph_p,dim=1).ne.mymdim).or.(size(graph_p,dim=2).ne.nats))then
2565  write(*,*)"PRG_COLLECT_GRAPH_P: graph_p has wrong shape. mdim=",mymdim,",nats=",nats,",shape=",shape(graph_p)
2566  stop
2567  endif
2568  endif
2569 
2570  lx = latticevectors(1,1)
2571  ly = latticevectors(2,2)
2572  lz = latticevectors(3,3)
2573 #ifdef USE_OFFLOAD
2574  allocate(graphed(nats,nc))
2575  allocate(rhoext(nats,nc))
2576  allocate(graph_core(mymdim,nc))
2577 
2578  rho_bml_c_ptr = bml_get_data_ptr_dense(rho_bml)
2579  ld = bml_get_ld_dense(rho_bml)
2580  call c_f_pointer(rho_bml_c_ptr,rho_bml_ptr,shape=[ld,norbs])
2581  !$acc enter data create(extmat(1:nats,1:nch),rho_red(1:nch,1:nc),graph_core(1:mymdim,1:nc)) &
2582  !$acc create(rhoext(1:nats,1:nc),graphed(1:nats,1:nc)) &
2583  !$acc copyin(hindex(1:2,1:nch),chindex(1:chindex_size))
2584 
2585  !$acc parallel loop gang &
2586  !$acc present(extmat)
2587 
2588  do i=1,nch
2589  !$acc loop vector private(dvx,dvy,dvz)
2590  do j = 1,nats
2591  dvx = modulo((coordsall(1,j) - coords(1,i) + lx/2.0_dp),lx) - lx/2.0_dp
2592  dvy = modulo((coordsall(2,j) - coords(2,i) + ly/2.0_dp),ly) - ly/2.0_dp
2593  dvz = modulo((coordsall(3,j) - coords(3,i) + lz/2.0_dp),lz) - lz/2.0_dp
2594  extmat(j,i) = exp(-alpha*(dvx*dvx+dvy*dvy+dvz*dvz))
2595  enddo
2596  !$acc end loop
2597  enddo
2598  !$acc end parallel loop
2599 
2600  !$acc parallel loop collapse(2) deviceptr(rho_bml_ptr) present(rho_red)
2601  do j=1,nc
2602  do i=1,nch
2603  rho_red(i,j) = maxval(abs(rho_bml_ptr(hindex(1,j):hindex(2,j),hindex(1,i):hindex(2,i))))
2604  enddo
2605  enddo
2606  !$acc end parallel loop
2607 
2608  !$acc parallel loop present(extmat,rho_red,rhoext) private(i,j,k)
2609  do i=1,nc
2610  do j=1,nats
2611  rhoext(j,i) = 0.0_dp
2612  do k=1,nch
2613  rhoext(j,i) = rhoext(j,i) + extmat(j,k)*rho_red(k,i)
2614  enddo
2615  enddo
2616  enddo
2617  !$acc end parallel loop
2618 
2619  !$acc parallel loop gang present(graphed,graph_core,rhoext) &
2620  !$acc present(rho_red,chindex) &
2621  !$acc private(ncounti,ii,i,ifull,j,jfull)
2622  do i = 1, nc
2623  ifull = chindex(i) + 1 !Map it to the full system
2624  !$acc loop vector
2625  do j = 1,nats
2626  graphed(j,i) = .false.
2627  enddo
2628  !$acc end loop
2629  !$acc loop
2630  do ii = 1,mymdim
2631  graph_core(ii,i) = 0
2632  enddo
2633  !$acc end loop
2634  !$acc loop
2635  do j = 1,nch
2636  if (rho_red(j,i).ge.threshold)then
2637  graphed(chindex(j)+1,i) = .true.
2638  endif
2639  enddo
2640  !$acc end loop
2641  ncounti = 0
2642  do j = 1, nats
2643  if ((rhoext(j,i).ge.threshold).or.graphed(j,i))then
2644  ncounti = ncounti + 1
2645  graph_core(ncounti,i) = j
2646  endif
2647  enddo
2648  enddo
2649  !$acc end parallel loop
2650 
2651  !$acc exit data copyout(graph_core(1:mymdim,1:nc)) &
2652  !$acc delete(graphed(1:nats,1:nc),rhoext(1:nats,1:nc),extmat(1:nats,1:nch)) &
2653  !$acc delete(rho_red(1:nch,1:nc),hindex(1:2,1:nch),chindex(1:chindex_size))
2654 
2655  !$omp parallel do shared(graph_p,graph_core,nc)
2656  do i = 1,nc
2657  graph_p(:,chindex(i) + 1) = graph_core(:,i) !Map it to the full system
2658  enddo
2659  !$omp end parallel do
2660 
2661  deallocate(graphed)
2662  deallocate(graph_core)
2663 #else
2664  allocate(rho(norbs,norbs))
2665 
2666  call bml_export_to_dense(rho_bml,rho)
2667 
2668  rho = abs(rho)
2669 
2670  !$omp parallel do default(none) &
2671  !$omp private(i,j) &
2672  !$omp private(dvec,dr2) &
2673  !$omp shared(extmat,alpha,coordsall,coords,Lx,Ly,Lz,nch,nc)
2674  do i=1,nch
2675  dvec(1,:) = modulo((coordsall(1,:) - coords(1,i) + lx/2.0_dp),lx) - lx/2.0_dp
2676  dvec(2,:) = modulo((coordsall(2,:) - coords(2,i) + ly/2.0_dp),ly) - ly/2.0_dp
2677  dvec(3,:) = modulo((coordsall(3,:) - coords(3,i) + lz/2.0_dp),lz) - lz/2.0_dp
2678  dr2(:) = dvec(1,:)*dvec(1,:) + dvec(2,:)*dvec(2,:) + dvec(3,:)*dvec(3,:)
2679  extmat(:,i) = exp(-alpha*dr2(:))
2680  enddo
2681  !$omp end parallel do
2682 
2683  !$omp parallel do collapse(2) default(none) &
2684  !$omp private(i,j) &
2685  !$omp shared(rho_red,rho,hindex,nch,nc)
2686  do j=1,nc
2687  do i=1,nch
2688  rho_red(i,j) = maxval(rho(hindex(1,i):hindex(2,i),hindex(1,j):hindex(2,j)))
2689  enddo
2690  enddo
2691  !$omp end parallel do
2692  allocate(rhoext(nats,nc))
2693  !rhoext = matmul(extmat,rho_red)
2694  !$omp parallel do private(i,j,k) shared(rhoext,extmat,rho_red) schedule(static)
2695  do i=1,nc
2696  do j=1,nats
2697  rhoext(j,i) = 0.0_dp
2698  do k=1,nch
2699  rhoext(j,i) = rhoext(j,i) + extmat(j,k)*rho_red(k,i)
2700  enddo
2701  enddo
2702  enddo
2703  !$omp end parallel do
2704 
2705  !$omp parallel do default(none) private(i) &
2706  !$omp private(ncounti,rowatfull,ii,j,jfull,ifull,iconnectedtoj) &
2707  !$omp private(row,connection,nats) &
2708  !$omp shared(graph_p,nc,chindex,hindex,rhoext,rho_red,nch,threshold)
2709  do i = 1, nc
2710  ifull = chindex(i) + 1 !Map it to the full system
2711  rowatfull = .false.
2712  ii=1
2713  ncounti = 0
2714  do while (graph_p(ii,ifull).ne.0) !Unfolding the connections of ifull
2715  rowatfull(graph_p(ii,ifull)) = .true.
2716  ncounti = ncounti + 1
2717  ii = ii+1
2718  enddo
2719  do j = 1,nch
2720  jfull = chindex(j) + 1 !Map it to the full system
2721  if ((rho_red(j,i).ge.threshold).and.(.not.rowatfull(jfull)))then
2722  ncounti = ncounti + 1
2723  graph_p(ncounti,ifull) = jfull
2724  rowatfull(jfull) = .true.
2725  endif
2726  enddo
2727  do j = 1, nats
2728  if ((rhoext(j,i).ge.threshold).and.(.not.rowatfull(j)))then
2729  ncounti = ncounti + 1
2730  graph_p(ncounti,ifull) = j
2731  endif
2732  enddo
2733  enddo
2734  !$omp end parallel do
2735  deallocate(rho)
2736 #endif
2737 
2738  deallocate(rowatfull)
2739  deallocate(iconnectedtoj)
2740  deallocate(row)
2741  deallocate(dvec)
2742  deallocate(dr2)
2743  deallocate(extmat)
2744  deallocate(rhoext)
2745 
2746  end subroutine prg_collect_extended_graph_p
2747 
2753  subroutine prg_merge_graph(graph_p,graph_h)
2754  implicit none
2755  integer :: i, ii, j, nats
2756  integer :: ncounti, maxnz
2757  integer, intent(inout) :: graph_h(:,:)
2758  integer, intent(inout) :: graph_p(:,:)
2759  logical, allocatable :: rowpatfull(:)
2760 
2761  nats = size(graph_p,dim=2)
2762  maxnz = size(graph_p,dim=1)
2763 
2764  if(.not.allocated(rowpatfull))allocate(rowpatfull(nats))
2765 
2766  !!$omp parallel do default(none) private(i) &
2767  !!$omp private(ncounti,rowpatfull,ii,j) &
2768  !!$omp shared(graph_p,graph_h,nats,maxnz)
2769  do i = 1, nats
2770  ncounti = 0
2771  rowpatfull = .false.
2772 
2773  ii=1
2774  do while ((graph_p(ii,i).ne.0).and.(ii.le.maxnz))
2775  if(graph_p(ii,i).gt.nats) then
2776  write(*,*)"PRG_MERGE_GRAPH: graph_p(",ii,",",i,") has value ",graph_p(ii,i)," which exceeds nats = ",nats
2777  stop
2778  endif
2779  rowpatfull(graph_p(ii,i)) = .true.
2780  ncounti = ncounti + 1
2781  ii = ii + 1
2782  enddo
2783 
2784  ii = 1
2785  do while ((graph_h(ii,i).ne.0).and.(ii.le.maxnz))
2786  j = graph_h(ii,i)
2787  if(.not.rowpatfull(j))then
2788  ncounti = ncounti + 1
2789  graph_p(ncounti,i) = j
2790  endif
2791  ii = ii + 1
2792  enddo
2793  enddo
2794  !!$omp end parallel do
2795 
2796  deallocate(rowpatfull)
2797 
2798  end subroutine prg_merge_graph
2799 
2800 
2808  subroutine prg_merge_graph_adj(graph_p,graph_h,xadj,adjncy)
2809  implicit none
2810  integer :: i, ii, j, nats
2811  integer :: ncounti, ncountot
2812  integer, allocatable, intent(inout) :: adjncy(:), graph_h(:,:), graph_p(:,:), xadj(:)
2813  logical(1), allocatable :: rowpatfull(:)
2814 
2815  nats = size(graph_p,dim=2)
2816  allocate(rowpatfull(nats))
2817  allocate(xadj(nats+1))
2818  allocate(adjncy(nats*nats))
2819 
2820  xadj = 0.0
2821  adjncy = 0.0
2822 
2823  !$omp parallel do default(none) private(i) &
2824  !$omp private(ncounti,rowpatfull,ii,j,ncountot) &
2825  !$omp shared(graph_p,graph_h,nats)
2826  do i = 1, nats
2827  ncounti = 0
2828  rowpatfull = .false.
2829  ii = 1
2830  do while (graph_p(ii,i).ne.0) !Unfolding the connections of i
2831  rowpatfull(graph_p(ii,i)) = .true.
2832  ncounti = ncounti + 1
2833  ii = ii+1
2834  enddo
2835 
2836  ii = 1
2837  do while (graph_h(ii,i).ne.0) !Unfolding the connections of i
2838  j = graph_h(ii,i)
2839  if(.not.rowpatfull(j))then
2840  ncounti = ncounti + 1
2841  ncountot = ncountot + 1
2842  graph_p(ncounti,i) = j
2843  endif
2844  ii = ii+1
2845  enddo
2846  enddo
2847  !$omp end parallel do
2848 
2849 
2850  ncountot = 0
2851 
2852  !$omp parallel do default(none) private(i) &
2853  !$omp private(ncounti,ii,ncountot) &
2854  !$omp shared(graph_p,nats,xadj,adjncy)
2855  do i = 1, nats
2856  ncounti = 0
2857  ii = 1
2858  do while (graph_p(ii,i).gt.0)
2859  ncounti = ncounti + 1
2860  ncountot = ncountot + 1
2861  adjncy(ncountot) = graph_p(ii,i)
2862  ii = ii+1
2863  enddo
2864  xadj(i) = ncountot - ncounti + 1
2865  enddo
2866  !$omp end parallel do
2867 
2868  xadj(nats+1) = ncountot + 1
2869 
2870  deallocate(rowpatfull)
2871  !deallocate(graph_p)
2872 
2873  end subroutine prg_merge_graph_adj
2874 
2882  subroutine prg_adj2bml(xadj,adjncy,bml_type,g_bml)
2883  implicit none
2884  character(20), intent(in) :: bml_type
2885  integer :: i, ii, j, nats
2886  integer, intent(in) :: adjncy(:), xadj(:)
2887  real(dp), allocatable :: row(:)
2888  type(bml_matrix_t), intent(inout) :: g_bml
2889 
2890  nats = size(xadj,dim=1) - 1
2891  if(bml_get_n(g_bml).gt.0) call bml_deallocate(g_bml)
2892  call bml_zero_matrix(bml_type,bml_element_real,dp,nats,nats,g_bml)
2893  allocate(row(nats))
2894 
2895  !$omp parallel do default(none) private(i) &
2896  !$omp private(j,row) &
2897  !$omp shared(nats,g_bml,xadj,adjncy)
2898  do i = 1, nats
2899  row = 0.0_dp
2900  do j = xadj(i), xadj(i+1) - 1
2901  row(adjncy(j)) = 1.0_dp
2902  enddo
2903  call bml_set_row(g_bml,i,row)
2904  enddo
2905 
2906  deallocate(row)
2907 
2908  end subroutine prg_adj2bml
2909 
2916  subroutine prg_graph2bml(graph,bml_type,g_bml)
2917  implicit none
2918  character(20), intent(in) :: bml_type
2919  integer :: i, ii, nats, mymdim
2920  integer, allocatable, intent(inout) :: graph(:,:)
2921  !real(dp), allocatable :: row(:)
2922  real(kind(1.0)), allocatable :: row(:)
2923  type(bml_matrix_t), intent(inout) :: g_bml
2924 
2925  nats = size(graph,dim=2)
2926  ! if(bml_get_N(g_bml).GT.0) call bml_deallocate(g_bml)
2927  ! call bml_zero_matrix(bml_type,bml_element_real,kind(1.0),nats,nats,g_bml)
2928  allocate(row(nats))
2929 
2930  mymdim = bml_get_m(g_bml)
2931  write(*,*)"mdim",mymdim
2932  !$omp parallel do default(none) private(i) &
2933  !$omp private(ii,row) &
2934  !$omp shared(nats,g_bml,graph,mymdim)
2935  do i = 1, nats
2936  ii = 1
2937  row = 0.0
2938  do while (graph(ii,i).gt.0)
2939  if(graph(ii,i).gt.nats)then
2940  write(*,*)"prg_graph2bml ERROR: nats < graph(",ii,",",i,") = ",graph(ii,i)
2941  stop
2942  endif
2943  row(graph(ii,i)) = 1.0_dp
2944  ! call bml_set_element_new(g_bml,i,graph(ii,i),1.0)
2945  ! call bml_set_element_new(g_bml,graph(ii,i),i,1.0)
2946  ii = ii+1
2947  enddo
2948  call bml_set_row(g_bml,i,row,0.5_dp)
2949  enddo
2950  !$omp end parallel do
2951 
2952  !deallocate(graph)
2953  deallocate(row)
2954 
2955  end subroutine prg_graph2bml
2956 
2957 
2963  subroutine prg_graph2vector(graph,vector,maxnz)
2964  implicit none
2965  integer, intent(inout) :: graph(:,:)
2966  integer, allocatable :: vector(:)
2967  integer :: i, j, maxnz, nats, ncount
2968 
2969  nats = size(graph,dim=2)
2970  allocate(vector(nats*maxnz))
2971 
2972  ncount = 0
2973 
2974  !$omp parallel do default(none) private(i) &
2975  !$omp private(ncount,j) &
2976  !$omp shared(nats,maxnz,vector,graph)
2977  do i = 1, nats
2978  do j = 1, maxnz
2979  ncount = (i-1)*maxnz + j !ncount + 1
2980  vector(ncount) = graph(j,i)
2981  enddo
2982  enddo
2983  !$omp end parallel do
2984 
2985  end subroutine prg_graph2vector
2986 
2992  subroutine prg_vector2graph(vector,graph,maxnz)
2993  implicit none
2994  integer, intent(inout) :: graph(:,:)
2995  integer, allocatable, intent(inout) :: vector(:)
2996  integer :: i, j, maxnz, nats, ncount
2997 
2998  nats = size(graph,dim=2)
2999  ncount = 0
3000 
3001  !$omp parallel do default(none) private(i) &
3002  !$omp private(ncount,j) &
3003  !$omp shared(nats,maxnz,vector,graph)
3004  do i = 1, nats
3005  do j = 1, maxnz
3006  ncount = (i-1)*maxnz + j !ncount + 1
3007  graph(j,i) = vector(ncount)
3008  enddo
3009  enddo
3010  !$omp end parallel do
3011 
3012  deallocate(vector)
3013 
3014  end subroutine prg_vector2graph
3015 
3020  subroutine prg_sortadj(xadj, adjncy)
3021  implicit none
3022  integer, intent(inout) :: xadj(:)
3023  integer, allocatable, intent(inout) :: adjncy(:)
3024  integer :: i, j, n, l, first, last, nx, irank, mintmp
3025  integer, allocatable :: tmpvect(:), newadjncy(:)
3026 
3027  nx = size(xadj,dim=1)
3028  allocate(tmpvect(nx))
3029  allocate(newadjncy(xadj(nx) - 1))
3030 
3031  !$omp parallel do default(none) private(i) &
3032  !$omp private(first,last,j,l,tmpvect,N,mintmp) &
3033  !$omp shared(nx,adjncy,xadj,newadjncy)
3034  do i = 1, nx-1
3035  first = xadj(i)
3036  last = xadj(i+1) - 1
3037  write(*,*)first,last
3038  n = last - first + 1
3039  tmpvect(1:n) = 0
3040  tmpvect(1:n) = adjncy(first:last)
3041 
3042  do j = 1,n
3043  do l = j+1,n
3044  if(tmpvect(j) >= tmpvect(l).and.tmpvect(l).ne.0)then
3045  mintmp = tmpvect(l)
3046  tmpvect(l) = tmpvect(j)
3047  tmpvect(j) = mintmp
3048  endif
3049  enddo
3050  enddo
3051 
3052  adjncy(first:last) = tmpvect(1:n)
3053  newadjncy(first:last) = tmpvect(1:n)
3054 
3055  enddo
3056  !$omp end parallel do
3057  deallocate(tmpvect)
3058  deallocate(adjncy)
3059 
3060  allocate(adjncy(xadj(nx) - 1)) !Reallocation for adjusting size
3061 
3062  adjncy = newadjncy
3063 
3064  deallocate(newadjncy)
3065 
3066  end subroutine prg_sortadj
3067 
3068 end module prg_system_mod
prg_extras_mod
Extra routines.
Definition: prg_extras_mod.F90:6
prg_extras_mod::dp
integer, parameter dp
Definition: prg_extras_mod.F90:31
prg_extras_mod::prg_norm2
real(dp) function, public prg_norm2(a)
Gets the norm2 of a vector.
Definition: prg_extras_mod.F90:241
prg_graph_mod
The graph module.
Definition: prg_graph_mod.F90:6
prg_graph_mod::prg_initgraphpartitioning
subroutine, public prg_initgraphpartitioning(gp, pname, np, nnodes, nnodes2)
Initialize graph partitioning.
Definition: prg_graph_mod.F90:181
prg_graph_mod::prg_initsubgraph
subroutine, public prg_initsubgraph(sg, pnum, hsize)
Initialize subgraph.
Definition: prg_graph_mod.F90:150
prg_graph_mod::prg_destroygraphpartitioning
subroutine, public prg_destroygraphpartitioning(gp)
Destroy graph partitioning.
Definition: prg_graph_mod.F90:269
prg_nvtx_mod
Definition: prg_nvtx_mod.F90:1
prg_openfiles_mod
Module to handle input output files for the PROGRESS lib.
Definition: prg_openfiles_mod.F90:4
prg_openfiles_mod::get_file_unit
integer function, public get_file_unit(io_max)
Returns a unit number that is not in use.
Definition: prg_openfiles_mod.F90:19
prg_openfiles_mod::prg_open_file
subroutine, public prg_open_file(io, name)
Opens a file to write.
Definition: prg_openfiles_mod.F90:38
prg_openfiles_mod::prg_open_file_to_read
subroutine, public prg_open_file_to_read(io, name)
Opens a file to read.
Definition: prg_openfiles_mod.F90:54
prg_ptable_mod
Periodic table of elements.
Definition: prg_ptable_mod.F90:8
prg_ptable_mod::element_covr
real(dp), dimension(nz), parameter element_covr
Covalent radius (in Angstroms)
Definition: prg_ptable_mod.F90:173
prg_ptable_mod::element_symbol
character(2), dimension(nz), parameter element_symbol
Element symbol.
Definition: prg_ptable_mod.F90:17
prg_ptable_mod::element_atomic_number
integer function, public element_atomic_number(symbol)
Definition: prg_ptable_mod.F90:394
prg_ptable_mod::element_atomic_number_upper
integer function element_atomic_number_upper(symbol)
Definition: prg_ptable_mod.F90:408
prg_ptable_mod::nz
integer, parameter nz
Definition: prg_ptable_mod.F90:12
prg_ptable_mod::element_mass
real(dp), dimension(nz), parameter element_mass
Element mass in atomic mass units (1.66 x 10-27 kg)
Definition: prg_ptable_mod.F90:110
prg_system_mod
A module to read and handle chemical systems.
Definition: prg_system_mod.F90:5
prg_system_mod::prg_get_subsystem
subroutine, public prg_get_subsystem(sy, lsize, indices, sbsy, verbose)
Get a subsystem out of the total system.
Definition: prg_system_mod.F90:2098
prg_system_mod::prg_adj2bml
subroutine, public prg_adj2bml(xadj, adjncy, bml_type, g_bml)
prg_adj2bml
Definition: prg_system_mod.F90:2883
prg_system_mod::prg_destroy_subsystems
subroutine, public prg_destroy_subsystems(sbsy, verbose)
Destroy allocated subsystem.
Definition: prg_system_mod.F90:2186
prg_system_mod::prg_parameters_to_vectors
subroutine, public prg_parameters_to_vectors(abc_angles, lattice_vector)
Transforms the lattice parameters into lattice vectors.
Definition: prg_system_mod.F90:1237
prg_system_mod::prg_translatetogeomcandfoldtobox
subroutine, public prg_translatetogeomcandfoldtobox(coords, lattice_vectors, origin)
Translate to geometric center.
Definition: prg_system_mod.F90:1514
prg_system_mod::prg_collect_extended_graph_p
subroutine, public prg_collect_extended_graph_p(rho_bml, nc, nats, hindex, chindex, graph_p, threshold, mdimin, alpha, coords, coordsall, latticevectors, verbose)
Collect the small graph to build the full graph.
Definition: prg_system_mod.F90:2503
prg_system_mod::prg_replicate_system
subroutine, public prg_replicate_system(sy, syf, nx, ny, nz)
Extend/replicate a system type along the lattice vectors.
Definition: prg_system_mod.F90:1610
prg_system_mod::prg_get_dihedral
subroutine, public prg_get_dihedral(coords, id1, id2, id3, id4, dihedral)
Get the dihedral angle given four atomic positions.
Definition: prg_system_mod.F90:1837
prg_system_mod::prg_graph2bml
subroutine, public prg_graph2bml(graph, bml_type, g_bml)
Graph2bml.
Definition: prg_system_mod.F90:2917
prg_system_mod::prg_write_trajectory
subroutine, public prg_write_trajectory(system, iter, each, prg_deltat, filename, extension)
Write trajectory in .xyz, .dat or pdb file.
Definition: prg_system_mod.F90:957
prg_system_mod::prg_get_covgraph_int
subroutine prg_get_covgraph_int(sy, nnStructMindist, nnStruct, nrnnstruct, bml_type, factor, gcov_bml, mdimin, verbose)
Definition: prg_system_mod.F90:1969
prg_system_mod::prg_molpartition
subroutine, public prg_molpartition(sy, npart, nnStructMindist, nnStruct, nrnnstruct, hetatm, gp, verbose)
Partition by molecule.
Definition: prg_system_mod.F90:2250
prg_system_mod::prg_make_random_system
subroutine, public prg_make_random_system(system, nats, seed, lx, ly, lz)
Make random Xx system.
Definition: prg_system_mod.F90:1191
prg_system_mod::prg_cleanuprepeatedatoms
subroutine, public prg_cleanuprepeatedatoms(nats, coords, symbols, verbose)
Cleanup repeated atoms we might have in the system.
Definition: prg_system_mod.F90:1707
prg_system_mod::prg_merge_graph_adj
subroutine, public prg_merge_graph_adj(graph_p, graph_h, xadj, adjncy)
Get partial subgraph based on the Density matrix.
Definition: prg_system_mod.F90:2809
prg_system_mod::prg_graph2vector
subroutine, public prg_graph2vector(graph, vector, maxnz)
Vectorize graph.
Definition: prg_system_mod.F90:2964
prg_system_mod::prg_get_covgraph
subroutine, public prg_get_covgraph(sy, nnStruct, nrnnstruct, bml_type, factor, gcov_bml, mdimin, verbose)
Get the covalency graph in bml format.
Definition: prg_system_mod.F90:1892
prg_system_mod::prg_write_trajectoryandproperty
subroutine, public prg_write_trajectoryandproperty(system, iter, each, prg_deltat, scalarprop, filename, extension)
Write trajectory and atomic properties. Only pdb file.
Definition: prg_system_mod.F90:1081
prg_system_mod::prg_vector2graph
subroutine, public prg_vector2graph(vector, graph, maxnz)
Back to graph.
Definition: prg_system_mod.F90:2993
prg_system_mod::prg_translateandfoldtobox
subroutine, public prg_translateandfoldtobox(coords, lattice_vectors, origin, verbose)
Translate and fold to box.
Definition: prg_system_mod.F90:1373
prg_system_mod::prg_wraparound
subroutine, public prg_wraparound(coords, lattice_vectors, index, verbose)
Wrap around atom i using pbc.
Definition: prg_system_mod.F90:1469
prg_system_mod::prg_destroy_estr
subroutine, public prg_destroy_estr(estr)
Deallocates all the arrays within the electronic structure.
Definition: prg_system_mod.F90:712
prg_system_mod::prg_get_covgraph_h
subroutine, public prg_get_covgraph_h(sy, nnStruct, nrnnstruct, rcut, graph_h, mdimin, verbose)
Get the covanlency graph.
Definition: prg_system_mod.F90:2029
prg_system_mod::prg_get_distancematrix
subroutine, public prg_get_distancematrix(coords, dmat)
Get the distance matrix.
Definition: prg_system_mod.F90:1348
prg_system_mod::prg_get_nameandext
subroutine, public prg_get_nameandext(fullfilename, filename, ext)
Get the name and extension of a file.
Definition: prg_system_mod.F90:251
prg_system_mod::prg_get_partial_atomgraph
subroutine, public prg_get_partial_atomgraph(rho_bml, hindex, gch_bml, threshold, verbose)
Get partial subgraph based on the Density matrix.
Definition: prg_system_mod.F90:2339
prg_system_mod::prg_destroy_system
subroutine, public prg_destroy_system(sy)
Deallocates all the arrays within a system.
Definition: prg_system_mod.F90:684
prg_system_mod::prg_parse_system
subroutine, public prg_parse_system(system, filename, extin)
The parser for the chemical system.
Definition: prg_system_mod.F90:274
prg_system_mod::prg_centeratbox
subroutine, public prg_centeratbox(coords, lattice_vectors, verbose)
Translate geometric center to the center of the box.
Definition: prg_system_mod.F90:1427
prg_system_mod::prg_vectors_to_parameters
subroutine, public prg_vectors_to_parameters(lattice_vector, abc_angles)
Transforms the lattice vectors into lattice parameters.
Definition: prg_system_mod.F90:1279
prg_system_mod::prg_merge_graph
subroutine, public prg_merge_graph(graph_p, graph_h)
Get partial subgraph based on the Density matrix.
Definition: prg_system_mod.F90:2754
prg_system_mod::prg_replicate
subroutine, public prg_replicate(coords, symbols, lattice_vectors, nx, ny, nz)
Extend/replicate system along lattice vectors.
Definition: prg_system_mod.F90:1553
prg_system_mod::prg_write_system
subroutine, public prg_write_system(system, filename, extin)
Write system in .xyz, .dat or pdb file.
Definition: prg_system_mod.F90:745
prg_system_mod::prg_get_recip_vects
subroutine, public prg_get_recip_vects(lattice_vectors, recip_vectors, volr, volk)
Get the volume of the cell and the reciprocal vectors: This soubroutine computes:
Definition: prg_system_mod.F90:1788
prg_system_mod::prg_sortadj
subroutine, public prg_sortadj(xadj, adjncy)
Sort adj NOTE: this might not be needed anymre since the bml_get_adj routine is sorting the values.
Definition: prg_system_mod.F90:3021
prg_system_mod::prg_get_origin
subroutine, public prg_get_origin(coords, origin)
Get the origin of the coordinates.
Definition: prg_system_mod.F90:1315
prg_system_mod::prg_collect_graph_p
subroutine, public prg_collect_graph_p(rho_bml, nc, nats, hindex, chindex, graph_p, threshold, mdimin, verbose)
Collect the small graph to build the full graph.
Definition: prg_system_mod.F90:2405
prg_graph_mod::graph_partitioning_t
Trace per iteration.
Definition: prg_graph_mod.F90:59
prg_system_mod::estruct_type
Electronic structure type.
Definition: prg_system_mod.F90:25
prg_system_mod::system_type
System type.
Definition: prg_system_mod.F90:112