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 prg_openfiles_mod
9  use prg_ptable_mod
10  use prg_graph_mod
11  use prg_extras_mod
12 
13  implicit none
14 
15  private
16 
17  integer, parameter :: dp = kind(1.0d0)
18 
20  type, public :: estruct_type
21 
23  integer :: norbs
24 
26  !approaches
27  integer :: norbscore
28 
30  integer :: nel
31 
33  integer, allocatable :: hindex(:,:)
34 
36  type(bml_matrix_t) :: ham
37 
39  type(bml_matrix_t) :: hamaux
40 
42  type(bml_matrix_t) :: ham0
43 
45  type(bml_matrix_t) :: oham
46 
48  type(bml_matrix_t) :: ohamaux
49 
51  type(bml_matrix_t) :: over
52 
54  type(bml_matrix_t) :: rho
55 
57  type(bml_matrix_t) :: orho
58 
60  type(bml_matrix_t) :: zmat
61 
63  type(bml_matrix_t) :: evects
64 
66  type(bml_matrix_t) :: mat
67 
69  type(bml_matrix_t) :: mataux
70 
72  real(dp), allocatable :: coul_pot_r(:)
73 
75  real(dp), allocatable :: coul_pot_k(:)
76 
78  real(dp), allocatable :: skforce(:,:)
79 
81  real(dp), allocatable :: fpul(:,:)
82 
84  real(dp), allocatable :: fscoul(:,:)
85 
87  real(dp), allocatable :: aux(:,:)
88 
90  real(dp), allocatable :: evals(:)
91 
93  real(dp), allocatable :: dvals(:)
94 
96  real(dp), allocatable :: ker(:,:)
97 
99  integer :: norbsincore
100 
102  real(dp) :: eband
103 
104  end type estruct_type
105 
107  type, public :: system_type
108 
110  integer :: nats = 0
111 
118  character(2), allocatable :: symbol(:)
119 
121  integer, allocatable :: atomic_number(:)
122 
126  real(dp), allocatable :: coordinate(:,:)
127 
131  real(dp), allocatable :: velocity(:,:)
132 
136  real(dp), allocatable :: force(:,:)
137 
141  real(dp), allocatable :: net_charge(:)
142 
148  real(dp), allocatable :: mass(:)
149 
158  real(dp), allocatable :: lattice_vector(:,:)
159 
166  real(dp), allocatable :: recip_vector(:,:)
167 
170  real(dp) :: volr
171 
174  real(dp) :: volk
175 
181  integer :: nsp
182 
189  integer, allocatable :: spindex(:)
190 
197  character(2), allocatable :: splist(:)
198 
204  integer, allocatable :: spatnum(:)
205 
210  real(dp), allocatable :: spmass(:)
211 
213  real(dp), allocatable :: userdef(:)
214 
216  integer, allocatable :: resindex(:)
217 
219  character(3), allocatable :: resname(:)
220 
222  character(3), allocatable :: atomname(:)
223 
225  type(estruct_type) :: estr
226 
227  end type system_type
228 
229  public :: prg_parse_system, prg_get_nameandext, prg_make_random_system, prg_write_system, prg_write_trajectory
230  public :: prg_get_origin, prg_get_covgraph, prg_get_subsystem, prg_translateandfoldtobox, prg_molpartition, prg_get_partial_atomgraph
231  public :: prg_destroy_subsystems, prg_get_covgraph_h, prg_collect_graph_p, prg_merge_graph, prg_merge_graph_adj, prg_adj2bml, prg_graph2bml
232  public :: prg_graph2vector, prg_vector2graph, prg_sortadj, prg_get_recip_vects, prg_translatetogeomcandfoldtobox
233  public :: prg_write_trajectoryandproperty, prg_get_distancematrix, prg_parameters_to_vectors, prg_vectors_to_parameters
234  public :: prg_get_dihedral, prg_wraparound, prg_centeratbox, prg_replicate, prg_cleanuprepeatedatoms
235  public :: prg_replicate_system, prg_destroy_system, prg_destroy_estr
236 
237 contains
238 
244  subroutine prg_get_nameandext(fullfilename,filename,ext)
245  implicit none
246  character(len=*), intent(in) :: fullfilename
247  character(50), intent(inout) :: filename
248  character(3), intent(inout) :: ext
249  character(1), allocatable :: tempc(:)
250  character(len=30) :: tempcflex
251  integer :: lenc
252 
253  lenc=len(adjustl(trim(fullfilename)))
254  if(.not.allocated(tempc))allocate(tempc(lenc))
255  tempcflex = adjustl(trim(fullfilename))
256  filename = adjustl(trim(tempcflex(1:lenc-4)))
257  ext = adjustl(trim(tempcflex(lenc-2:lenc+1)))
258 
259  end subroutine prg_get_nameandext
260 
261 
267  subroutine prg_parse_system(system,filename,extin)
268  implicit none
269  character(1) :: onechar
270  character(10) :: dummyc(10)
271  character(2) :: possiblenewsymbol, twochar
272  character(2), allocatable :: sptempsymbols(:)
273  character(3), optional, intent(in) :: extin
274  character(3) :: extension
275  character(30) :: dummy
276  character(50) :: io_name, nametmp
277  character(60) :: pdbformat
278  character(len=*) :: filename
279  integer :: dummyi(10), header_lines, i, io_unit
280  integer :: ios, itemp, j, lines_to_lattice
281  integer :: max_lines, nats, nsp, verbose
282  integer :: lines_to_masses, lines_to_atom
283  logical :: islattice
284  real(dp) :: abc_angles(2,3), dummyr(10), lbx1
285  real(dp) :: lbx2, lby1, lby2, lbz1
286  real(dp) :: lbz2, max_x, max_y, max_z
287  real(dp) :: min_x, min_y, min_z, scfactor
288  type(system_type), intent(out) :: system
289 
290  verbose = 0
291  islattice = .false.
292 
293  max_x = 0.0_dp; max_y = 0.0_dp; max_z = 0.0_dp;
294  min_x = 0.0_dp; min_y = 0.0_dp; min_z = 0.0_dp;
295 
296  if(allocated(system%symbol)) stop "ERROR: System already allocated"
297 
298  if(.not.present(extin))then
299  call prg_get_nameandext(filename,nametmp,extension)
300  else
301  extension = extin
302  nametmp = trim(adjustl(filename))
303  endif
304 
305  select case(extension)
306 
307  case("xyz")
308 
309  !! For xyz format see http://openbabel.org/wiki/XYZ_%28format%29
310  io_name=trim(adjustl(nametmp))//".xyz"
311  call prg_open_file_to_read(io_unit,io_name)
312  read(io_unit,*)nats
313  read(io_unit,*)
314  system%nats = nats
315  allocate(system%symbol(nats))
316  allocate(system%atomic_number(nats))
317  allocate(system%coordinate(3,nats))
318  allocate(system%mass(nats))
319  allocate(system%lattice_vector(3,3))
320 
321  system%lattice_vector = 0.0_dp
322 
323  do i=1,nats
324  read(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i)
325  system%atomic_number(i) = element_atomic_number(system%symbol(i))
326  system%mass(i) = element_mass(system%atomic_number(i))
327  max_x = max(max_x,system%coordinate(1,i))
328  min_x = min(min_x,system%coordinate(1,i))
329  max_y = max(max_y,system%coordinate(2,i))
330  min_y = min(min_y,system%coordinate(2,i))
331  max_z = max(max_z,system%coordinate(3,i))
332  min_z = min(min_z,system%coordinate(3,i))
333  ! write(*,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i)
334  enddo
335 
336  !The following is not part of an xyz format.
337  !VMD, babel, xmakemol and pymol can still read this file.
338  read(io_unit,*,iostat=ios)dummy
339  if(dummy.eq."#lattice")then
340  if(verbose.eq.1)write(*,*)"There is a lattice ..."
341  read(io_unit,*)system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
342  read(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
343  read(io_unit,*)system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
344  else
345  !We will create the lattice vectors if they are not pressent.
346  !It will add 10.0 Ang to each coordinate.
347  system%lattice_vector = 0.0_dp
348  system%lattice_vector(1,1) = max_x - min_x + 10.0_dp
349  system%lattice_vector(2,2) = max_y - min_y + 10.0_dp
350  system%lattice_vector(3,3) = max_z - min_z + 10.0_dp
351  endif
352 
353  close(io_unit)
354 
355  case("pdb")
356 
357  !! For PDB format see http://www.wwpdb.org/documentation/file-format
358  io_name=trim(adjustl(nametmp))//".pdb"
359  call prg_open_file_to_read(io_unit,io_name)
360  header_lines = 0
361  lines_to_lattice = 0
362  max_lines = 1000000
363 
364  !Counting header lines
365  do i=1,max_lines
366  read(io_unit,*)dummy
367  if(dummy.eq."ATOM".or.dummy.eq."HETATM") then
368  exit
369  else
370  header_lines = header_lines + 1
371  if(dummy.eq."CRYST1")then
372  lines_to_lattice = header_lines
373  islattice = .true.
374  endif
375  endif
376  enddo
377 
378  nats = 1
379  do i=1,max_lines
380  read(io_unit,*)dummy
381  if(dummy.eq."ATOM".or.dummy.eq."HETATM") then
382  nats = nats + 1
383  else
384  exit
385  endif
386  enddo
387  close(io_unit)
388 
389  call prg_open_file_to_read(io_unit,io_name)
390  do i=1,header_lines
391  if(i.eq.lines_to_lattice)then
392  read(io_unit,*)dummy,abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
393  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)
394  else
395  read(io_unit,*)dummy
396  endif
397  enddo
398 
399  system%nats = nats
400  allocate(system%symbol(nats))
401  allocate(system%atomic_number(nats))
402  allocate(system%coordinate(3,nats))
403  allocate(system%mass(nats))
404  allocate(system%lattice_vector(3,3))
405  allocate(system%resname(nats))
406  allocate(system%resindex(nats))
407  allocate(system%atomname(nats))
408 
409  system%lattice_vector = 0.0_dp
410 
411  pdbformat= '(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2)'
412 
413  do i=1,nats
414  read(io_unit,pdbformat)dummyc(1),dummyi(1), &
415  system%atomname(i),dummyc(3),system%resname(i),dummyc(4),system%resindex(i),dummyc(5),&
416  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
417  dummyr(1),dummyr(2),system%symbol(i),dummyc(10)
418 
419  ! In case there are no symbols in the last column:
420  if(dummyc(4).ne."".and.system%symbol(i).eq."")then
421  onechar=adjustl(trim(dummyc(4)))
422  if(onechar.ne."H".and. &
423  onechar.ne."B".and. &
424  onechar.ne."C".and. &
425  onechar.ne."N".and. &
426  onechar.ne."O".and. &
427  onechar.ne."F".and. &
428  onechar.ne."P".and. &
429  onechar.ne."S".and. &
430  onechar.ne."I")then
431  twochar=adjustl(trim(dummyc(4)))
432  system%symbol(i)=twochar
433  else
434  system%symbol(i)=onechar
435  endif
436  system%atomic_number(i) = element_atomic_number_upper(system%symbol(i))
437  system%symbol(i) = element_symbol(system%atomic_number(i)) !upper to lower char.
438  else
439  system%atomic_number(i) = element_atomic_number(system%symbol(i))
440  endif
441 
442  ! Getting the atomic mass
443  system%mass(i) = element_mass(system%atomic_number(i))
444 
445  ! Getting the system limits.
446  max_x = max(max_x,system%coordinate(1,i))
447  min_x = min(min_x,system%coordinate(1,i))
448  max_y = max(max_y,system%coordinate(2,i))
449  min_y = min(min_y,system%coordinate(2,i))
450  max_z = max(max_z,system%coordinate(3,i))
451  min_z = min(min_z,system%coordinate(3,i))
452 
453  enddo
454 
455  if(islattice)then
456  call prg_parameters_to_vectors(abc_angles,system%lattice_vector)
457  else
458  !We will create the lattice vectors if they are not pressent
459  !It will add 10.0 Ang to each coordinate.
460  system%lattice_vector = 0.0_dp
461  system%lattice_vector(1,1) = max_x - min_x + 10.0_dp
462  system%lattice_vector(2,2) = max_y - min_y + 10.0_dp
463  system%lattice_vector(3,3) = max_z - min_z + 10.0_dp
464  endif
465 
466  close(io_unit)
467 
468  case("ltt")
469 
470  !! For old inputblock.dat (old dat) format see the LATTE manual.
471  io_name=trim(adjustl(nametmp))//".ltt"
472  call prg_open_file_to_read(io_unit,io_name)
473  read(io_unit,*)dummy, nats
474  read(io_unit,*)scfactor
475  system%nats = nats
476  allocate(system%symbol(nats))
477  allocate(system%atomic_number(nats))
478  allocate(system%coordinate(3,nats))
479  allocate(system%mass(nats))
480  allocate(system%lattice_vector(3,3))
481 
482  !!The gets the llattice vectors from the box boundaries.
483  read(io_unit,*)lbx1,lbx2,lby1,lby2,lbz1,lbz2
484  system%lattice_vector = 0.0_dp
485  system%lattice_vector(1,1) = (lbx2-lbx1)*scfactor
486  system%lattice_vector(2,2) = (lby2-lby1)*scfactor
487  system%lattice_vector(3,3) = (lbz2-lbz1)*scfactor
488 
489  do i=1,nats
490  read(io_unit,*)system%symbol(i),system%coordinate(1,i)&
491  ,system%coordinate(2,i),system%coordinate(3,i)
492  system%coordinate(1,i)= system%coordinate(1,i)*scfactor
493  system%coordinate(2,i)= system%coordinate(2,i)*scfactor
494  system%coordinate(3,i)= system%coordinate(3,i)*scfactor
495  system%atomic_number(i) = element_atomic_number(system%symbol(i))
496  system%mass(i) = element_mass(system%atomic_number(i))
497  enddo
498 
499  close(io_unit)
500 
501  case("dat")
502 
503  !! For new inputblock.dat (dat) format see the LATTE manual.
504  io_name=trim(adjustl(nametmp))//".dat"
505  call prg_open_file_to_read(io_unit,io_name)
506  read(io_unit,*)nats
507  system%nats = nats
508  allocate(system%symbol(nats))
509  allocate(system%atomic_number(nats))
510  allocate(system%coordinate(3,nats))
511  allocate(system%mass(nats))
512  allocate(system%lattice_vector(3,3))
513 
514  !!The gets the llattice vectors from the box boundaries.
515  system%lattice_vector = 0.0_dp
516  read(io_unit,*)system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
517  read(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
518  read(io_unit,*)system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
519 
520  do i=1,nats
521  read(io_unit,*)system%symbol(i),system%coordinate(1,i)&
522  ,system%coordinate(2,i),system%coordinate(3,i)
523  system%coordinate(1,i)= system%coordinate(1,i)
524  system%coordinate(2,i)= system%coordinate(2,i)
525  system%coordinate(3,i)= system%coordinate(3,i)
526  system%atomic_number(i) = element_atomic_number(system%symbol(i))
527  system%mass(i) = element_mass(system%atomic_number(i))
528  enddo
529 
530  close(io_unit)
531 
532  case("gen")
533 
534  stop "gen format is not implemented for reading yet"
535 
536  case("lmp")
537 
538  !! For Lammps data.* input file.
539  !! For more information see: http://lammps.sandia.gov/doc/2001/data_format.html
540  io_name=adjustl(trim(nametmp))//".lmp"
541  call prg_open_file(io_unit,io_name)
542  read(io_unit,*)dummyc(1)
543  read(io_unit,*)system%nats
544 
545  do i=1,100
546  read(io_unit,*)dummyc(1),dummyc(2)!,dummyc(3)
547 
548  if(adjustl(trim(dummyc(2))) == "atom")then
549  lines_to_atom = i + 2
550  endif
551  if(adjustl(trim(dummyc(1))) == "Masses")then
552  lines_to_masses = i + 2
553  exit
554  endif
555  enddo
556  close(io_unit)
557 
558  call prg_open_file(io_unit,io_name)
559 
560  do i=1,lines_to_atom-1
561  read(io_unit,*)dummyc(1)
562  enddo
563 
564  read(io_unit,*)system%nsp
565 
566  do i=1,lines_to_masses - lines_to_atom - 4
567  read(io_unit,*)dummyc(1)
568  enddo
569 
570  read(io_unit,*)min_x, max_x
571  read(io_unit,*)min_y, max_y
572  read(io_unit,*)min_z, max_z
573 
574  allocate(system%lattice_vector(3,3))
575  system%lattice_vector = 0.0_dp
576  system%lattice_vector(1,1) = max_x - min_x
577  system%lattice_vector(2,2) = max_y - min_y
578  system%lattice_vector(3,3) = max_z - min_z
579 
580  read(io_unit,*)dummyc(1)
581 
582  allocate(system%spmass(system%nsp))
583  allocate(system%splist(system%nsp))
584 
585  do i=1,system%nsp
586  read(io_unit,*)dummyi(1),system%spmass(i)
587  enddo
588 
589  do i=1,system%nsp
590  do j=1,size(element_mass,dim=1)
591  if(abs(system%spmass(i) - element_mass(j)) < 1.0) then
592  system%splist(i) = element_symbol(j)
593  exit
594  endif
595  enddo
596  enddo
597 
598  do i=1,100000
599  read(io_unit,*)dummyc(1)
600  if(dummyc(1) == "Atoms")exit
601  enddo
602 
603  allocate(system%spindex(system%nats))
604  allocate(system%coordinate(3,system%nats))
605  allocate(system%symbol(system%nats))
606 
607  do i=1,system%nats
608  read(io_unit,*)dummyi(1),dummyi(2),system%spindex(i),dummyr(1),system%coordinate(1,i)&
609  ,system%coordinate(2,i),system%coordinate(3,i)
610  write(*,*)dummyi(1),dummyi(2),system%spindex(i),dummyr(1),system%coordinate(1,i)&
611  ,system%coordinate(2,i),system%coordinate(3,i)
612  system%symbol(i) = system%splist(system%spindex(i))
613  enddo
614 
615  close(io_unit)
616 
617  case default
618 
619  stop "The file extension is not valid. Only xyz, lmp, dat and pdb formats are implemented"
620 
621  end select
622 
623  !!Computing species and species index
624  !!Assigning species index to the system elements.
625  nsp=1
626  allocate(sptempsymbols(150)) !Temporary vector to accumulate the species.
627  sptempsymbols(1)=system%Symbol(1)
628 
629  if(.not.allocated(system%spindex)) allocate(system%spindex(nats))
630 
631  do i=1,system%nats
632  itemp = 0
633  do j=1,nsp !If there is a new chemical symbol
634  if(trim(system%symbol(i)).ne.trim(sptempsymbols(j)))then
635  itemp = itemp + 1
636  possiblenewsymbol = system%symbol(i)
637  endif
638  enddo
639  if(itemp == nsp)then
640  nsp = nsp+1
641  sptempsymbols(nsp)=possiblenewsymbol
642  endif
643  enddo
644 
645  if(allocated(system%splist))deallocate(system%splist)
646  if(allocated(system%spmass))deallocate(system%spmass)
647 
648  allocate(system%splist(nsp))
649  allocate(system%spmass(nsp))
650  allocate(system%spatnum(nsp))
651 
652  do i=1,nsp
653  system%splist(i) = sptempsymbols(i)
654  system%spatnum(i) = element_atomic_number(system%splist(i))
655  system%spmass(i) = element_mass(system%spatnum(i))
656  enddo
657 
658  deallocate(sptempsymbols)
659  system%nsp = nsp
660 
663  do i = 1,system%nats
664  do j=1,nsp
665  if(trim(system%symbol(i)).eq.trim(system%splist(j)))then
666  system%spindex(i)=j
667  endif
668  enddo
669  enddo
670 
671  end subroutine prg_parse_system
672 
673 
677  subroutine prg_destroy_system(sy)
678  implicit none
679  type(system_type), intent(inout) :: sy
680 
681  if(allocated(sy%symbol)) deallocate(sy%symbol)
682  if(allocated(sy%atomic_number)) deallocate(sy%atomic_number)
683  if(allocated(sy%coordinate)) deallocate(sy%coordinate)
684  if(allocated(sy%velocity)) deallocate(sy%velocity)
685  if(allocated(sy%force)) deallocate(sy%force)
686  if(allocated(sy%net_charge)) deallocate(sy%net_charge)
687  if(allocated(sy%mass)) deallocate(sy%mass)
688  if(allocated(sy%lattice_vector)) deallocate(sy%lattice_vector)
689  if(allocated(sy%recip_vector)) deallocate(sy%recip_vector)
690  if(allocated(sy%spindex)) deallocate(sy%spindex)
691  if(allocated(sy%splist)) deallocate(sy%splist)
692  if(allocated(sy%spatnum)) deallocate(sy%spatnum)
693  if(allocated(sy%spmass)) deallocate(sy%spmass)
694  if(allocated(sy%userdef)) deallocate(sy%userdef)
695  if(allocated(sy%resindex)) deallocate(sy%resindex)
696  if(allocated(sy%resname)) deallocate(sy%resname)
697  if(allocated(sy%atomname)) deallocate(sy%atomname)
698 
699  end subroutine prg_destroy_system
700 
704  subroutine prg_destroy_estr(estr)
705  implicit none
706  type(estruct_type), intent(inout) :: estr
707 
708  if(allocated(estr%hindex))deallocate(estr%hindex)
709  if(bml_allocated(estr%ham))call bml_deallocate(estr%ham)
710  if(bml_allocated(estr%hamAux)) call bml_deallocate(estr%hamAux)
711  if(bml_allocated(estr%ham0)) call bml_deallocate(estr%ham0)
712  if(bml_allocated(estr%oham)) call bml_deallocate(estr%oham)
713  if(bml_allocated(estr%ohamAux)) call bml_deallocate(estr%ohamAux)
714  if(bml_allocated(estr%over)) call bml_deallocate(estr%over)
715  if(bml_allocated(estr%rho)) call bml_deallocate(estr%rho)
716  if(bml_allocated(estr%evects)) call bml_deallocate(estr%evects)
717  if(bml_allocated(estr%matAux)) call bml_deallocate(estr%matAux)
718  if(bml_allocated(estr%zmat)) call bml_deallocate(estr%zmat)
719  if(allocated(estr%coul_pot_r)) deallocate(estr%coul_pot_r)
720  if(allocated(estr%coul_pot_k)) deallocate(estr%coul_pot_k)
721  if(allocated(estr%skforce)) deallocate(estr%skforce)
722  if(allocated(estr%fpul)) deallocate(estr%fpul)
723  if(allocated(estr%fscoul)) deallocate(estr%fscoul)
724  if(allocated(estr%aux)) deallocate(estr%aux)
725  if(allocated(estr%evals)) deallocate(estr%evals)
726  if(allocated(estr%dvals)) deallocate(estr%dvals)
727  if(allocated(estr%ker)) deallocate(estr%ker)
728 
729  end subroutine prg_destroy_estr
730 
731 
737  subroutine prg_write_system(system,filename,extin)
738  implicit none
739  character(len=*) :: filename
740  character(10) :: dummyc(10)
741  character(34) :: xyzformat
742  character(3), optional, intent(in) :: extin
743  character(3) :: extension
744  character(50) :: io_name, nametmp
745  character(60) :: pdbformat
746  character(100) :: io_message
747  integer :: dummyi(10), i, io_unit, nats
748  real(dp) :: abc_angles(2,3), dummyr(10), origin(3)
749  type(system_type), intent(inout) :: system
750 
751  nats = system%nats
752 
753  if(.not.present(extin))then
754  call prg_get_nameandext(filename,nametmp,extension)
755  else
756  extension = extin
757  nametmp = trim(filename)
758  endif
759 
760  select case(extension)
761 
762  case("xyz")
763 
764  io_name=trim(nametmp)//".xyz"
765  call prg_open_file(io_unit,io_name)
766  write(io_unit,*)nats
767  io_message = trim(adjustl(io_name))//" Generated by the PROGRESS library"
768  write(io_unit,*)trim(adjustl(io_message))
769  xyzformat = '(A2,A2,1F10.5,A2,1F10.5,A2,1F10.5)'
770  do i=1,nats
771  write(io_unit,xyzformat)system%symbol(i)," ",system%coordinate(1,i)," ",system%coordinate(2,i)," ",system%coordinate(3,i)
772  enddo
773 
774  !The following is not part of an xyz format but
775  !VMD, babel, xmakemol and pymol can still read this file.
776  write(io_unit,*)"#lattice vectors"
777  write(io_unit,"(3F20.5)")system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
778  write(io_unit,"(3F20.5)")system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
779  write(io_unit,"(3F20.5)")system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
780 
781  close(io_unit)
782 
783  case("pdb")
784 
785  dummyi = 0
786  dummyc = ""
787  dummyc(1) = "ATOM"
788  dummyr = 0.0_dp
789 
790  io_name=trim(nametmp)//".pdb"
791  call prg_open_file(io_unit,io_name)
792 
793  write(io_unit,'(A6,1X,A50)')"REMARK","Generated by PROGRESS library"
794  write(io_unit,'(A5,1X,A20)')"TITLE",io_name
795 
796  if(allocated(system%lattice_vector))then
797  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
798  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
799  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
800  endif
801 
802  write(io_unit,'(A5,A20)')"MODEL"," 1"
803 
804  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2 )'
805 
806  if(.not.allocated(system%resindex))then
807  allocate(system%resindex(nats))
808  system%resindex = 1
809  endif
810 
811  if(.not.allocated(system%resname))then
812  allocate(system%resname(nats))
813  system%resname = "MOL"
814  endif
815 
816  if(.not.allocated(system%atomname))then
817  allocate(system%atomname(nats))
818  system%atomname = system%symbol
819  endif
820 
821  if(allocated(system%net_charge))then
822  do i=1,nats
823  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
824  system%atomname(i),dummyc(5),system%resname(i),dummyc(7),system%resindex(i),dummyc(8),&
825  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
826  dummyr(1),system%net_charge(i),system%symbol(i),dummyc(10)
827  enddo
828  else
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),dummyr(2),system%symbol(i),dummyc(10)
834  enddo
835  endif
836 
837  write(io_unit,'(A3)')"TER"
838  write(io_unit,'(A3)')"ENDMDL"
839  close(io_unit)
840 
841 
842  case("ltt")
843 
844  !! For old inputblock.dat (old dat) format see the LATTE manual.
845  io_name=trim(nametmp)//".ltt"
846  call prg_open_file(io_unit,io_name)
847  write(io_unit,*)"NATS= ", system%nats
848 
849  !Here we will always write the system in 1:1 scale
850  write(io_unit,*)"1.0"," Generated by the PROGRESS library"
851 
852  !! This gets the lattice vectors from the box boundaries.
853  write(io_unit,"(6F10.5)")0.0,system%lattice_vector(1,1),0.0,system%lattice_vector(2,2)&
854  ,0.0,system%lattice_vector(3,3)
855 
856  do i=1,nats
857  write(io_unit,"(A2,3F10.5)")system%symbol(i),system%coordinate(1,i)&
858  ,system%coordinate(2,i),system%coordinate(3,i)
859  enddo
860 
861  close(io_unit)
862 
863  case("dat")
864 
865  !! For new inputblock.dat (nat) format see the LATTE manual.
866  io_name=trim(nametmp)//".dat"
867  call prg_open_file(io_unit,io_name)
868  write(io_unit,*)system%nats
869 
870  !! This gets the lattice vectors from the box boundaries.
871  write(io_unit,"(3F10.5)")system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
872  write(io_unit,"(3F10.5)")system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
873  write(io_unit,"(3F10.5)")system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
874 
875  do i=1,nats
876  write(io_unit,"(A2,3F10.5)")system%symbol(i),system%coordinate(1,i)&
877  ,system%coordinate(2,i),system%coordinate(3,i)
878  enddo
879 
880  close(io_unit)
881 
882  case("gen")
883 
884  !! For gen format see DFTB+ manual.
885  io_name=trim(nametmp)//".gen"
886  call prg_open_file(io_unit,io_name)
887  write(io_unit,*)system%nats,"S"
888  write(io_unit,'(103A3)')(system%splist(i),i=1,system%nsp)
889  do i=1,nats
890  write(io_unit,"(I10,I10,3F10.5)")i,system%spindex(i),system%coordinate(1,i)&
891  ,system%coordinate(2,i),system%coordinate(3,i)
892  enddo
893  write(io_unit,*)"0.0 0.0 0.0"
894  write(io_unit,'(3F10.5)')system%lattice_vector(1,1),system%lattice_vector(1,2),system%lattice_vector(1,3)
895  write(io_unit,'(3F10.5)')system%lattice_vector(2,1),system%lattice_vector(2,2),system%lattice_vector(2,3)
896  write(io_unit,'(3F10.5)')system%lattice_vector(3,1),system%lattice_vector(3,2),system%lattice_vector(3,3)
897 
898  close(io_unit)
899 
900  case("lmp")
901 
902  !! For Lammps data.* input file.
903  !! For more information see: http://lammps.sandia.gov/doc/2001/data_format.html
904  io_name= adjustl(trim(nametmp))//".lmp"
905  call prg_open_file(io_unit,io_name)
906  write(io_unit,*)"LAMMPS Description"
907  write(io_unit,*)""
908  write(io_unit,*)system%nats,"atoms"
909  write(io_unit,*)""
910  write(io_unit,*)system%nsp,"atom types"
911  write(io_unit,*)""
912  write(io_unit,*)0.0_dp,system%lattice_vector(1,1),"xlo xhi"
913  write(io_unit,*)0.0_dp,system%lattice_vector(2,2),"ylo yhi"
914  write(io_unit,*)0.0_dp,system%lattice_vector(3,3),"zlo zhi"
915  write(io_unit,*)system%lattice_vector(2,1),system%lattice_vector(3,1), &
916  &system%lattice_vector(3,2), "xy xz yz"
917 
918  write(io_unit,*)""
919  write(io_unit,*)"Masses"
920  write(io_unit,*)""
921  do i=1,system%nsp
922  write(io_unit,*)" ",i,system%spmass(i)
923  enddo
924  write(io_unit,*)""
925  write(io_unit,*)"Atoms"
926  write(io_unit,*)""
927  do i=1,nats
928  write(io_unit,'(3I5,A6,3F10.5)')i,1,system%spindex(i)," 0.0",system%coordinate(1,i)&
929  ,system%coordinate(2,i),system%coordinate(3,i)
930  enddo
931 
932  close(io_unit)
933 
934  case default
935 
936  stop "The file extension is not valid. Only xyz, dat, lmp, pdb or gen formats are implemented"
937 
938  end select
939 
940  end subroutine prg_write_system
941 
949  subroutine prg_write_trajectory(system,iter,each,prg_deltat,filename,extension)
950  implicit none
951  character(len=*) :: filename
952  character(10) :: dummyc(10)
953  character(11) :: xyzformat
954  character(20) :: io_name
955  character(3) :: extension
956  character(60) :: pdbformat
957  integer :: dummyi(10), i, io_unit, nats
958  integer, intent(in) :: iter, each
959  integer, allocatable :: resindex(:)
960  real(dp) :: abc_angles(2,3), dummyr(10)
961  real(dp), intent(in) :: prg_deltat
962  type(system_type), intent(in) :: system
963 
964  if(mod(iter,each).eq.0.or.iter.eq.0.or.iter.eq.1)then
965  nats = system%nats
966 
967  select case(extension)
968 
969  case("xyz")
970 
971  io_name=trim(filename)//".xyz"
972  io_unit=get_file_unit(-1)
973 
974  if(iter.eq.0.or.iter.eq.1)then
975  open(unit=io_unit,file=io_name,status='unknown')
976  else
977  open(unit=io_unit,file=io_name,position = 'append',status='old')
978  endif
979 
980  write(io_unit,*)system%nats
981  write(io_unit,*)"frame", iter
982  if(allocated(system%net_charge))then
983  do i=1,nats
984  write(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
985  system%net_charge(i)
986  enddo
987  else
988  do i=1,nats
989  write(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i)
990  enddo
991  endif
992 
993  close(io_unit)
994 
995  case("pdb")
996 
997  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
998 
999  dummyi = 0
1000  dummyc = ""
1001  dummyc(1) = "ATOM"
1002  dummyr = 0.0_dp
1003 
1004  io_name=trim(filename)//".pdb"
1005  io_unit=get_file_unit(-1)
1006 
1007  if(iter.eq.0.or.iter.eq.1)then
1008  open(unit=io_unit,file=io_name,status='unknown')
1009  else
1010  open(unit=io_unit,file=io_name,position = 'append',status='old')
1011  endif
1012 
1013  write(io_unit,'(A17,A4,F10.5)')"REMARK Trajectory"," t=",iter*prg_deltat
1014  write(io_unit,'(A31)')"REMARK THIS IS A SIMULATION BOX"
1015  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
1016  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
1017  write(io_unit,'(A5,I6)')"MODEL",iter
1018 
1019  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2 )'
1020 
1021  if(.not.allocated(system%resindex))then
1022  allocate(resindex(nats))
1023  resindex = 1
1024  else
1025  resindex = system%resindex
1026  endif
1027 
1028 
1029  if(allocated(system%net_charge))then
1030  do i=1,nats
1031  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1032  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1033  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1034  dummyr(1),system%net_charge(i),system%symbol(i),dummyc(10)
1035  enddo
1036  else
1037  do i=1,nats
1038  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1039  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1040  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1041  dummyr(1),dummyr(2),system%symbol(i),dummyc(10)
1042  enddo
1043  endif
1044 
1045  write(io_unit,'(A3)')"TER"
1046  write(io_unit,'(A6)')"ENDMDL"
1047  close(io_unit)
1048 
1049  case("dat")
1050 
1051  write(*,*)"write traj not implemented for dat "
1052 
1053  case default
1054 
1055  stop "The file extension is not valid. Only pdb and xyz formats are implemented"
1056 
1057  end select
1058 
1059  endif
1060 
1061  end subroutine prg_write_trajectory
1062 
1072  subroutine prg_write_trajectoryandproperty(system,iter,each,prg_deltat,scalarprop,filename,extension)
1073  implicit none
1074  character(len=*) :: filename
1075  character(10) :: dummyc(10)
1076  character(11) :: xyzformat
1077  character(20) :: io_name
1078  character(3) :: extension
1079  character(60) :: pdbformat
1080  integer :: dummyi(10), i, io_unit, nats
1081  integer, intent(in) :: iter, each
1082  integer, allocatable :: resindex(:)
1083  real(dp) :: abc_angles(2,3), dummyr(10)
1084  real(dp) :: maxprop, minprop, realtmp
1085  real(dp), intent(in) :: prg_deltat, scalarprop(:)
1086  type(system_type), intent(in) :: system
1087 
1088 
1089  ! maxprop = maxval(scalarprop(:))
1090  ! minprop = minval(scalarprop(:))
1091 
1092  if(mod(iter,each).eq.0.or.iter.eq.0.or.iter.eq.1)then
1093  nats = system%nats
1094 
1095  select case(extension)
1096 
1097  case("pdb")
1098 
1099  call prg_vectors_to_parameters(system%lattice_vector,abc_angles)
1100 
1101  dummyi = 0
1102  dummyc = ""
1103  dummyc(1) = "ATOM"
1104  dummyr = 0.0_dp
1105 
1106  io_name=trim(filename)//".pdb"
1107  io_unit=get_file_unit(-1)
1108 
1109  if(iter.eq.0.or.iter.eq.1)then
1110  open(unit=io_unit,file=io_name,status='unknown')
1111  else
1112  open(unit=io_unit,file=io_name,position = 'append',status='old')
1113  endif
1114 
1115  write(io_unit,'(A8,A4,F10.5)')"Trajectory"," t=",iter*prg_deltat
1116  write(io_unit,'(A24)')"THIS IS A SIMULATION BOX"
1117  write(io_unit,'(A6,3F9.3,3F7.2,A16)')"CRYST1",abc_angles(1,1),abc_angles(1,2),abc_angles(1,3)&
1118  ,abc_angles(2,1),abc_angles(2,2),abc_angles(2,3)," P 1 1"
1119  write(io_unit,'(A5,I6)')"MODEL",iter
1120 
1121  pdbformat= '(A4,A2,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2 )'
1122 
1123  if(.not.allocated(system%resindex))then
1124  allocate(resindex(nats))
1125  resindex = 1
1126  else
1127  resindex = system%resindex
1128  endif
1129 
1130  do i=1,nats
1131  ! realtmp = (scalarprop(i) - minprop)/(maxprop-minprop) - 1.0_dp !Normalization
1132  write(io_unit,pdbformat)dummyc(1),dummyc(2),i, &
1133  system%symbol(i),dummyc(5),"MOL",dummyc(7),resindex(i),dummyc(8),&
1134  system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1135  dummyr(1),scalarprop(i),system%symbol(i),dummyc(10)
1136  enddo
1137 
1138  write(io_unit,'(A3)')"TER"
1139  write(io_unit,'(A3)')"ENDMDL"
1140  close(io_unit)
1141 
1142  case("xyz")
1143 
1144  io_name=trim(filename)//".xyz"
1145  io_unit=get_file_unit(-1)
1146 
1147  if(iter.eq.0.or.iter.eq.1)then
1148  open(unit=io_unit,file=io_name,status='unknown')
1149  else
1150  open(unit=io_unit,file=io_name,position = 'append',status='old')
1151  endif
1152 
1153  write(io_unit,*)system%nats
1154  write(io_unit,*)"frame", iter
1155  do i=1,nats
1156  ! realtmp = (scalarprop(i) - minprop)/(maxprop-minprop) - 1.0_dp !Normalization
1157  write(io_unit,*)system%symbol(i),system%coordinate(1,i),system%coordinate(2,i),system%coordinate(3,i),&
1158  scalarprop(i)
1159  ! write(*,*)scalarprop(i)
1160  enddo
1161 
1162  close(io_unit)
1163 
1164  case default
1165 
1166  stop "The file extension is not valid. Only pdb format is implemented"
1167 
1168  end select
1169 
1170  endif
1171 
1172  end subroutine prg_write_trajectoryandproperty
1173 
1174 
1182  subroutine prg_make_random_system(system,nats,seed,lx,ly,lz)
1183 
1184  implicit none
1185  integer :: i, nats, seed
1186  real(dp) :: lx, ly, lz, ran
1187  type(system_type), intent(out) :: system
1188 
1189  !Allocating the system
1190  system%nats = nats
1191  allocate(system%symbol(nats))
1192  allocate(system%atomic_number(nats))
1193  allocate(system%coordinate(3,nats))
1194  allocate(system%mass(nats))
1195 
1196  call random_seed(seed)
1197 
1198  do i=1,nats
1199 
1200  system%atomic_number(i) = 1
1201  call random_number(ran)
1202  system%coordinate(i,1) = ran*lx
1203  call random_number(ran)
1204  system%coordinate(i,2) = ran*ly
1205  call random_number(ran)
1206  system%coordinate(i,3) = ran*lz
1207 
1208  system%symbol(i) = element_symbol(system%atomic_number(i))
1209  system%mass(i) = element_mass(1)
1210 
1211  enddo
1212 
1213  allocate(system%lattice_vector(3,3))
1214  system%lattice_vector = 0.0_dp
1215  system%lattice_vector(1,1) = lx
1216  system%lattice_vector(2,2) = ly
1217  system%lattice_vector(3,3) = lz
1218 
1219  end subroutine prg_make_random_system
1220 
1228  subroutine prg_parameters_to_vectors(abc_angles,lattice_vector)
1229  implicit none
1230  real(dp) :: angle_alpha, angle_beta, angle_gamma, lattice_a
1231  real(dp) :: lattice_b, lattice_c, pi
1232  real(dp), intent(in) :: abc_angles(2,3)
1233  real(dp), intent(out) :: lattice_vector(3,3)
1234 
1235  pi = 3.14159265358979323846264338327950_dp
1236 
1237  lattice_a = abc_angles(1,1)
1238  lattice_b = abc_angles(1,2)
1239  lattice_c = abc_angles(1,3)
1240  angle_alpha = abc_angles(2,1)
1241  angle_beta = abc_angles(2,2)
1242  angle_gamma = abc_angles(2,3)
1243 
1244  angle_alpha = 2.0_dp*pi*angle_alpha/360.0_dp
1245  angle_beta = 2.0_dp*pi*angle_beta/360.0_dp
1246  angle_gamma = 2.0_dp*pi*angle_gamma/360.0_dp
1247 
1248  lattice_vector(1,1) = lattice_a
1249  lattice_vector(1,2)=0
1250  lattice_vector(1,3)=0
1251 
1252  lattice_vector(2,1)=lattice_b*cos(angle_gamma)
1253  lattice_vector(2,2)=lattice_b*sin(angle_gamma)
1254  lattice_vector(2,3)=0
1255 
1256  lattice_vector(3,1)=lattice_c*cos(angle_beta)
1257  lattice_vector(3,2)=lattice_c*( cos(angle_alpha)-cos(angle_gamma)* &
1258  cos(angle_beta) )/sin(angle_gamma)
1259  lattice_vector(3,3)=sqrt(lattice_c**2 - lattice_vector(3,1)**2 - lattice_vector(3,2)**2)
1260 
1261  end subroutine prg_parameters_to_vectors
1262 
1270  subroutine prg_vectors_to_parameters(lattice_vector,abc_angles)
1271  implicit none
1272  real(dp) :: angle_alpha, angle_beta, angle_gamma, lattice_a
1273  real(dp) :: lattice_b, lattice_c, pi
1274  real(dp), intent(in) :: lattice_vector(3,3)
1275  real(dp), intent(out) :: abc_angles(2,3)
1276 
1277  pi = 3.14159265358979323846264338327950_dp
1278 
1279  lattice_a = sqrt(lattice_vector(1,1)**2+lattice_vector(1,2)**2+lattice_vector(1,3)**2)
1280  lattice_b = sqrt(lattice_vector(2,1)**2+lattice_vector(2,2)**2+lattice_vector(2,3)**2)
1281  lattice_c = sqrt(lattice_vector(3,1)**2+lattice_vector(3,2)**2+lattice_vector(3,3)**2)
1282 
1283  angle_gamma = dot_product(lattice_vector(1,:),lattice_vector(2,:))/(lattice_a*lattice_b)
1284  angle_beta = dot_product(lattice_vector(1,:),lattice_vector(3,:))/(lattice_a*lattice_c)
1285  angle_alpha = dot_product(lattice_vector(2,:),lattice_vector(3,:))/(lattice_b*lattice_c)
1286 
1287 
1288  angle_alpha = 360.0_dp*acos(angle_alpha)/(2.0_dp*pi)
1289  angle_beta = 360.0_dp*acos(angle_beta)/(2.0_dp*pi)
1290  angle_gamma = 360.0_dp*acos(angle_gamma)/(2.0_dp*pi)
1291 
1292  abc_angles(1,1) = lattice_a
1293  abc_angles(1,2) = lattice_b
1294  abc_angles(1,3) = lattice_c
1295 
1296  abc_angles(2,1) = angle_alpha
1297  abc_angles(2,2) = angle_beta
1298  abc_angles(2,3) = angle_gamma
1299 
1300  end subroutine prg_vectors_to_parameters
1301 
1306  subroutine prg_get_origin(coords,origin)
1307  implicit none
1308  integer :: i
1309  real(dp) :: max_x, max_y, max_z, min_x
1310  real(dp) :: min_y, min_z
1311  real(dp), allocatable, intent(inout) :: origin(:)
1312  real(dp), intent(in) :: coords(:,:)
1313 
1314  if(.not.allocated(origin)) allocate(origin(3))
1315 
1316  max_x = -1.0d5 ; max_y = -1.0d5 ; max_z = -1.0d5 ;
1317  min_x = 1.0d5 ; min_y = 1.0d5 ; min_z = 1.0d5 ;
1318 
1319  ! Getting the system limits.
1320  do i=1,size(coords,dim=2)
1321  max_x = max(max_x,coords(1,i))
1322  min_x = min(min_x,coords(1,i))
1323  max_y = max(max_y,coords(2,i))
1324  min_y = min(min_y,coords(2,i))
1325  max_z = max(max_z,coords(3,i))
1326  min_z = min(min_z,coords(3,i))
1327  enddo
1328 
1329  origin(1) = min_x
1330  origin(2) = min_y
1331  origin(3) = min_z
1332 
1333  end subroutine prg_get_origin
1334 
1339  subroutine prg_get_distancematrix(coords,dmat)
1340  implicit none
1341  integer :: i,j,nats
1342  real(dp) :: norm2
1343  real(dp), intent(in) :: coords(:,:)
1344  real(dp), intent(out), allocatable :: dmat(:,:)
1345 
1346  nats = size(coords,dim=2)
1347  if(.not.allocated(dmat)) allocate(dmat(nats,nats))
1348 
1349  ! Getting the system limits.
1350  do i=1,nats
1351  do j=1,nats
1352  dmat(i,j) = prg_norm2(coords(:,i)-coords(:,j))
1353  enddo
1354  enddo
1355 
1356  end subroutine prg_get_distancematrix
1357 
1358 
1364  subroutine prg_translateandfoldtobox(coords,lattice_vectors,origin, verbose)
1365  implicit none
1366  integer :: i
1367  integer, intent(in), optional :: verbose
1368  real(dp) :: max_x, max_y, max_z, min_x
1369  real(dp) :: min_y, min_z
1370  real(dp), allocatable, intent(inout) :: origin(:),coords(:,:)
1371  real(dp), intent(in) :: lattice_vectors(:,:)
1372 
1373  if(present(verbose))then
1374  if(verbose >= 1)write(*,*)"In prg_translateandfoldtobox ..."
1375  endif
1376 
1377  if(.not.allocated(origin)) allocate(origin(3))
1378 
1379  max_x = -1.0d5 ; max_y = -1.0d5 ; max_z = -1.0d5 ;
1380  min_x = 1.0d5 ; min_y = 1.0d5 ; min_z = 1.0d5 ;
1381 
1382  ! Getting the system limits.
1383  do i=1,size(coords,dim=2)
1384  max_x = max(max_x,coords(1,i))
1385  min_x = min(min_x,coords(1,i))
1386  max_y = max(max_y,coords(2,i))
1387  min_y = min(min_y,coords(2,i))
1388  max_z = max(max_z,coords(3,i))
1389  min_z = min(min_z,coords(3,i))
1390  enddo
1391 
1392  origin(1) = min_x
1393  origin(2) = min_y
1394  origin(3) = min_z
1395 
1396  coords(1,:) = coords(1,:) - origin(1)
1397  coords(2,:) = coords(2,:) - origin(2)
1398  coords(3,:) = coords(3,:) - origin(3)
1399 
1400  coords(1,:) = mod(coords(1,:)*1000000,1000000*lattice_vectors(1,1))/1000000;
1401  coords(2,:) = mod(coords(2,:)*1000000,1000000*lattice_vectors(2,2))/1000000;
1402  coords(3,:) = mod(coords(3,:)*1000000,1000000*lattice_vectors(3,3))/1000000;
1403 
1404  origin = 0.0_dp
1405 
1406  !Get new origin.
1407  origin(1) = -1.0d-1 ; origin(2) = -1.0d-1; origin(3) = -1.0d-1
1408 
1409  end subroutine prg_translateandfoldtobox
1410 
1411 
1417  subroutine prg_centeratbox(coords,lattice_vectors,verbose)
1418  implicit none
1419  integer :: i, nats
1420  integer, intent(in), optional :: verbose
1421  real(dp) :: gc(3)
1422  real(dp), allocatable, intent(inout) :: coords(:,:)
1423  real(dp), intent(in) :: lattice_vectors(:,:)
1424 
1425  if(present(verbose))then
1426  if(verbose >= 1)write(*,*)"In prg_centeratbox ..."
1427  endif
1428 
1429  nats=size(coords,dim=2)
1430 
1431  gc= 0.0d0
1432 
1433  !$omp parallel do default(none) private(i) &
1434  !$omp shared(coords,nats) &
1435  !$omp reduction(+:gc)
1436  do i=1,nats
1437  gc=gc + coords(:,i)
1438  enddo
1439  !$omp end parallel do
1440 
1441  gc=gc/real(nats,dp)
1442 
1443  !$omp parallel do default(none) private(i) &
1444  !$omp shared(coords,lattice_vectors,nats, gc)
1445  do i=1,nats
1446  coords(1,i) = coords(1,i) + lattice_vectors(1,1)/2.0d0 - gc(1)
1447  coords(2,i) = coords(2,i) + lattice_vectors(2,2)/2.0d0 - gc(2)
1448  coords(3,i) = coords(3,i) + lattice_vectors(3,3)/2.0d0 - gc(3)
1449  enddo
1450  !$omp end parallel do
1451 
1452  end subroutine prg_centeratbox
1453 
1459  subroutine prg_wraparound(coords,lattice_vectors,index,verbose)
1460  implicit none
1461  integer :: i, nats
1462  integer, intent(in) :: index
1463  integer, intent(in), optional :: verbose
1464  real(dp), allocatable, intent(inout) :: coords(:,:)
1465  real(dp), allocatable :: origin(:)
1466  real(dp), intent(in) :: lattice_vectors(:,:)
1467 
1468  if(present(verbose))then
1469  if(verbose >= 1)write(*,*)"In prg_wraparound ..."
1470  endif
1471 
1472  if(.not.allocated(origin)) allocate(origin(3))
1473 
1474  nats=size(coords,dim=2)
1475 
1476  origin(1) = -coords(1,index) + lattice_vectors(1,1)/2.0_dp
1477  origin(2) = -coords(2,index) + lattice_vectors(2,2)/2.0_dp
1478  origin(3) = -coords(3,index) + lattice_vectors(3,3)/2.0_dp
1479 
1480  coords(1,:) = coords(1,:) + origin(1)
1481  coords(2,:) = coords(2,:) + origin(2)
1482  coords(3,:) = coords(3,:) + origin(3)
1483 
1484  !$omp parallel do default(none) private(i) &
1485  !$omp shared(coords,lattice_vectors,nats)
1486  do i=1,nats
1487  if(coords(1,i) > lattice_vectors(1,1))coords(1,i)=coords(1,i)-lattice_vectors(1,1)
1488  if(coords(2,i) > lattice_vectors(2,2))coords(2,i)=coords(2,i)-lattice_vectors(2,2)
1489  if(coords(3,i) > lattice_vectors(3,3))coords(3,i)=coords(3,i)-lattice_vectors(3,3)
1490  if(coords(1,i) < 0.0_dp)coords(1,i)=coords(1,i)+lattice_vectors(1,1)
1491  if(coords(2,i) < 0.0_dp)coords(2,i)=coords(2,i)+lattice_vectors(2,2)
1492  if(coords(3,i) < 0.0_dp)coords(3,i)=coords(3,i)+lattice_vectors(3,3)
1493  enddo
1494  !$end omp parallel do
1495 
1496  end subroutine prg_wraparound
1497 
1498 
1504  subroutine prg_translatetogeomcandfoldtobox(coords,lattice_vectors,origin)
1505  implicit none
1506  integer :: i
1507  real(dp), allocatable, intent(inout) :: origin(:),coords(:,:)
1508  real(dp), intent(in) :: lattice_vectors(:,:)
1509  real(dp) :: geomc(3)
1510 
1511  if(.not.allocated(origin)) allocate(origin(3))
1512 
1513  ! Getting the geometric center.
1514  do i=1,size(coords,dim=2)
1515  geomc = geomc + coords(:,i)
1516  enddo
1517 
1518  geomc = geomc/size(coords,dim=2)
1519 
1520  coords(1,:) = coords(1,:) - geomc(1)
1521  coords(2,:) = coords(2,:) - geomc(2)
1522  coords(3,:) = coords(3,:) - geomc(3)
1523 
1524  coords(1,:) = mod(coords(1,:)*1000000,1000000*lattice_vectors(1,1))/1000000;
1525  coords(2,:) = mod(coords(2,:)*1000000,1000000*lattice_vectors(2,2))/1000000;
1526  coords(3,:) = mod(coords(3,:)*1000000,1000000*lattice_vectors(3,3))/1000000;
1527 
1528  origin = 0.0_dp
1529 
1530  !Shift origin slightly
1531  origin(1) = -1.0d-1 ; origin(2) = -1.0d-1; origin(3) = -1.0d-1
1532 
1533  end subroutine prg_translatetogeomcandfoldtobox
1534 
1543  subroutine prg_replicate(coords,symbols,lattice_vectors,nx,ny,nz)
1544  implicit none
1545  integer :: i, j, k, l
1546  integer :: m, fnats, nats
1547  integer, intent(in) :: nx, ny, nz
1548  character(2), allocatable, intent(inout) :: symbols(:)
1549  character(2), allocatable :: rsymbols(:)
1550  real(dp), allocatable, intent(inout) :: coords(:,:)
1551  real(dp), intent(inout) :: lattice_vectors(:,:)
1552  real(dp), allocatable :: r(:,:)
1553 
1554  nats = size(coords,dim=2)
1555  fnats = nx*ny*nz*nats
1556 
1557  allocate(r(3,fnats))
1558  allocate(rsymbols(fnats))
1559 
1560  m = 0
1561  do i=1,nx
1562  do j=1,ny
1563  do k=1,nz
1564  do l=1,nats
1565  m=m+1
1566  r(1,m)=i*lattice_vectors(1,1) + j*lattice_vectors(1,2) + k*lattice_vectors(1,3)
1567  r(1,m)=r(1,m) + coords(1,l)
1568  r(2,m)=i*lattice_vectors(2,1) + j*lattice_vectors(2,2) + k*lattice_vectors(2,3)
1569  r(2,m)=r(2,m) + coords(2,l)
1570  r(3,m)=i*lattice_vectors(3,1) + j*lattice_vectors(3,2) + k*lattice_vectors(3,3)
1571  r(3,m)=r(3,m) + coords(3,l)
1572  rsymbols(m) = symbols(l)
1573  enddo
1574  enddo
1575  enddo
1576  enddo
1577 
1578  deallocate(coords,symbols)
1579  allocate(coords(3,fnats))
1580  allocate(symbols(fnats))
1581 
1582  lattice_vectors(1,:) = nx*lattice_vectors(1,:)
1583  lattice_vectors(2,:) = ny*lattice_vectors(2,:)
1584  lattice_vectors(3,:) = nz*lattice_vectors(3,:)
1585 
1586  coords = r
1587  symbols = rsymbols
1588  deallocate(r,rsymbols)
1589 
1590  end subroutine prg_replicate
1591 
1592 
1600  subroutine prg_replicate_system(sy,syf,nx,ny,nz)
1601  implicit none
1602  integer :: i, j, k, l
1603  integer :: ini, fin, cont, ii
1604  integer, intent(in) :: nx, ny, nz
1605  type(system_type), intent(inout) :: sy
1606  type(system_type) :: syf
1607 
1608  if(.not. allocated(sy%coordinate))then
1609  write(*,*)"ERROR: System does not have coordinate"
1610  stop
1611  endif
1612  if(.not. allocated(sy%symbol))then
1613  write(*,*)"ERROR: System does not have atomic symbols"
1614  stop
1615  endif
1616  syf%nats = sy%nats*nx*ny*nz
1617  allocate(syf%coordinate(3,syf%nats))
1618  allocate(syf%symbol(syf%nats))
1619  allocate(syf%resindex(syf%nats))
1620  allocate(syf%atomname(syf%nats))
1621  allocate(syf%atomic_number(syf%nats))
1622  allocate(syf%spindex(syf%nats))
1623 
1624  cont = 0
1625  do i = 1,nx
1626  do j = 1,ny
1627  do k = 1,nz
1628  do ii = 1,sy%nats
1629  cont = cont + 1
1630  syf%coordinate(1,cont) = sy%coordinate(1,ii) + &
1631  & sy%lattice_vector(1,1)*i + sy%lattice_vector(2,1)*i + sy%lattice_vector(3,1)*i
1632  syf%coordinate(2,cont) = sy%coordinate(2,ii) + &
1633  & sy%lattice_vector(1,2)*j + sy%lattice_vector(2,2)*j + sy%lattice_vector(3,2)*j
1634  syf%coordinate(3,cont) = sy%coordinate(3,ii) + &
1635  & sy%lattice_vector(1,3)*k + sy%lattice_vector(2,3)*k + sy%lattice_vector(3,3)*k
1636  syf%symbol(cont) = sy%symbol(ii)
1637  enddo
1638  enddo
1639  enddo
1640  enddo
1641 
1642  if(allocated(sy%resindex))then
1643  do i = 1,nx*ny*nz
1644  ini = (i - 1)*sy%nats + 1
1645  fin = i*sy%nats
1646  syf%resindex(ini:fin) = sy%resindex(:)
1647  enddo
1648  else
1649  syf%resindex = 1
1650  endif
1651 
1652  if(allocated(sy%atomname))then
1653  do i = 1,nx*ny*nz
1654  ini = (i - 1)*sy%nats + 1
1655  fin = i*sy%nats
1656  syf%atomname(ini:fin) = sy%atomname(:)
1657  enddo
1658  else
1659  syf%atomname = "At"
1660  endif
1661 
1662  if(allocated(sy%atomic_number))then
1663  do i = 1,nx*ny*nz
1664  ini = (i - 1)*sy%nats + 1
1665  fin = i*sy%nats
1666  syf%atomic_number(ini:fin) = sy%atomic_number(:)
1667  enddo
1668  endif
1669  if(allocated(sy%spindex))then
1670  do i = 1,nx*ny*nz
1671  ini = (i - 1)*sy%nats + 1
1672  fin = i*sy%nats
1673  syf%spindex(ini:fin) = sy%spindex(:)
1674  enddo
1675  endif
1676 
1677  allocate(syf%lattice_vector(3,3))
1678  syf%lattice_vector(1,:) = nx*sy%lattice_vector(1,:)
1679  syf%lattice_vector(2,:) = ny*sy%lattice_vector(2,:)
1680  syf%lattice_vector(3,:) = nz*sy%lattice_vector(3,:)
1681 
1682  if(allocated(sy%splist))then
1683  syf%nsp = sy%nsp
1684  allocate(syf%splist(syf%nsp))
1685  syf%splist = sy%splist
1686  endif
1687 
1688  end subroutine prg_replicate_system
1689 
1690 
1697  subroutine prg_cleanuprepeatedatoms(nats,coords,symbols,verbose)
1698  implicit none
1699  integer :: i, natstmp, l, count, j, count1
1700  integer, intent(inout) :: nats
1701  integer, intent(in), optional :: verbose
1702  real(dp) :: d
1703  real(dp), allocatable, intent(inout) :: coords(:,:)
1704  real(dp), allocatable :: coordstmp(:,:)
1705  character(len=*), allocatable, intent(inout) :: symbols(:)
1706  character(2), allocatable :: symbolstmp(:)
1707 
1708  if(present(verbose))then
1709  if(verbose >= 1)write(*,*)"In prg_cleanuprepeatedatoms ..."
1710  endif
1711 
1712  natstmp=nats
1713 
1714  allocate(coordstmp(3,nats))
1715  allocate(symbolstmp(nats))
1716 
1717  coordstmp(1,1)=coords(1,1)
1718  coordstmp(2,1)=coords(2,1)
1719  coordstmp(3,1)=coords(3,1)
1720 
1721  l=1
1722  symbolstmp(1)=symbols(1)
1723 
1724  do j=1,nats
1725  count=0
1726  count1=0
1727  do i=1,l
1728  d=sqrt((coords(1,j)-coordstmp(1,i))**2+(coords(2,j)-coordstmp(2,i))**2+&
1729  &(coords(3,j)-coordstmp(3,i))**2)
1730  if(d.gt.0.1d0)then
1731  count=count + 1
1732  else
1733  count1=count1 + 1
1734  endif
1735  enddo
1736 
1737  if(count.eq.l)then ! New atom in collection
1738  if(count1 >= 1)then
1739  write(*,*)"Problem at iteration",j
1740  stop
1741  endif
1742 
1743  l = l+1
1744  coordstmp(1,l)=coords(1,j)
1745  coordstmp(2,l)=coords(2,j)
1746  coordstmp(3,l)=coords(3,j)
1747  symbolstmp(l)=symbols(j)
1748  endif
1749  enddo
1750 
1751  deallocate(symbols)
1752  deallocate(coords)
1753  nats = l
1754  allocate(symbols(nats))
1755  allocate(coords(3,nats))
1756 
1757  do i=1,nats
1758  symbols(i) = symbolstmp(i)
1759  coords(:,i) = coordstmp(:,i)
1760  enddo
1761 
1762  deallocate(coordstmp)
1763  deallocate(symbolstmp)
1764 
1765  end subroutine prg_cleanuprepeatedatoms
1766 
1778  subroutine prg_get_recip_vects(lattice_vectors,recip_vectors,volr,volk)
1779  implicit none
1780  real(dp) :: a1xa2(3), a2xa3(3), a3xa1(3), b2xb3(3)
1781  real(dp) :: pi
1782  real(dp), allocatable, intent(inout) :: recip_vectors(:,:)
1783  real(dp), intent(in) :: lattice_vectors(:,:)
1784  real(dp), intent(inout) :: volk, volr
1785 
1786  if(.not.allocated(recip_vectors))allocate(recip_vectors(3,3))
1787  volr=0.0_dp; volk=0.0_dp
1788 
1789  pi = 3.14159265358979323846264338327950_dp
1790 
1791  a1xa2(1) = lattice_vectors(1,2)*lattice_vectors(2,3) - lattice_vectors(1,3)*lattice_vectors(2,2)
1792  a1xa2(2) = -lattice_vectors(1,1)*lattice_vectors(2,3) + lattice_vectors(1,3)*lattice_vectors(2,1)
1793  a1xa2(3) = lattice_vectors(1,1)*lattice_vectors(2,2) - lattice_vectors(1,2)*lattice_vectors(2,1)
1794 
1795  a2xa3(1) = lattice_vectors(2,2)*lattice_vectors(3,3) - lattice_vectors(2,3)*lattice_vectors(3,2)
1796  a2xa3(2) = -lattice_vectors(2,1)*lattice_vectors(3,3) + lattice_vectors(2,3)*lattice_vectors(3,1)
1797  a2xa3(3) = lattice_vectors(2,1)*lattice_vectors(3,2) - lattice_vectors(2,2)*lattice_vectors(3,1)
1798 
1799  a3xa1(1) = lattice_vectors(3,2)*lattice_vectors(1,3) - lattice_vectors(3,3)*lattice_vectors(1,2)
1800  a3xa1(2) = -lattice_vectors(3,1)*lattice_vectors(1,3) + lattice_vectors(3,3)*lattice_vectors(1,1)
1801  a3xa1(3) = lattice_vectors(3,1)*lattice_vectors(1,2) - lattice_vectors(3,2)*lattice_vectors(1,1)
1802 
1803  !Get the volume of the cell
1804  volr = lattice_vectors(1,1)*a2xa3(1)+ lattice_vectors(1,2)*a2xa3(2)+lattice_vectors(1,3)*a2xa3(3)
1805 
1806  !Get the reciprocal vectors
1807  recip_vectors(1,:) = 2.0_dp*pi*a2xa3(:)/volr
1808  recip_vectors(2,:) = 2.0_dp*pi*a3xa1(:)/volr
1809  recip_vectors(3,:) = 2.0_dp*pi*a1xa2(:)/volr
1810 
1811  b2xb3(1) = recip_vectors(2,2)*recip_vectors(3,3) - recip_vectors(2,3)*recip_vectors(3,2)
1812  b2xb3(2) = -recip_vectors(2,1)*recip_vectors(3,3) + recip_vectors(2,3)*recip_vectors(3,1)
1813  b2xb3(3) = recip_vectors(2,1)*recip_vectors(3,2) - recip_vectors(2,2)*recip_vectors(3,1)
1814 
1815  volk = recip_vectors(1,1)*b2xb3(1)+ recip_vectors(1,2)*b2xb3(2)+recip_vectors(1,3)*b2xb3(3)
1816 
1817  end subroutine prg_get_recip_vects
1818 
1827  subroutine prg_get_dihedral(coords,id1,id2,id3,id4,dihedral)
1828 
1829  real(dp) :: mv1, mv2, v1(3), v2(3)
1830  real(dp) :: dotprod, cosdir, v2xv20(3), v1xv10(3)
1831  real(dp) :: v10(3),v20(3), cprod(3), normcprod, sindir
1832  real(dp), intent(in) :: coords(:,:)
1833  real(dp), intent(out) :: dihedral
1834  integer :: i
1835  integer, intent(in) :: id1,id2,id3,id4
1836  character(2) :: index1, index2, index3, index4
1837 
1838  v1=coords(:,id4) - coords(:,id3)
1839  v10=coords(:,id2) - coords(:,id3)
1840  v2=coords(:,id1) - coords(:,id2)
1841  v20=coords(:,id3) - coords(:,id2)
1842 
1843  v1xv10(1)=v1(2)*v10(3)-v1(3)*v10(2)
1844  v1xv10(2)=-(v1(1)*v10(3)-v1(3)*v10(1))
1845  v1xv10(3)=v1(1)*v10(2)-v1(2)*v10(1)
1846 
1847  v2xv20(1)=v2(2)*v20(3)-v2(3)*v20(2)
1848  v2xv20(2)=-(v2(1)*v20(3)-v2(3)*v20(1))
1849  v2xv20(3)=v2(1)*v20(2)-v2(2)*v20(1)
1850 
1851  dotprod = v1xv10(1)*v2xv20(1) + v1xv10(2)*v2xv20(2) + v1xv10(3)*v2xv20(3)
1852 
1853  cprod(1)=v1xv10(2)*v2xv20(3)-v1xv10(3)*v2xv20(2)
1854  cprod(2)=-(v1xv10(1)*v2xv20(3)-v1xv10(3)*v2xv20(1))
1855  cprod(3)=v1xv10(1)*v2xv20(2)-v1xv10(2)*v2xv20(1)
1856 
1857  normcprod=sqrt(cprod(1)*cprod(1) + cprod(2)*cprod(2) + cprod(3)*cprod(3))
1858 
1859  mv1= sqrt(v1xv10(1)*v1xv10(1) + v1xv10(2)*v1xv10(2) + v1xv10(3)*v1xv10(3))
1860  mv2= sqrt(v2xv20(1)*v2xv20(1) + v2xv20(2)*v2xv20(2) + v2xv20(3)*v2xv20(3))
1861 
1862  cosdir = dotprod/(mv1*mv2)
1863  sindir = normcprod/(mv1*mv2)
1864 
1865  sindir = cprod(3)
1866  dihedral=sign(1.0d0,sindir)*acos(-cosdir)
1867  dihedral=-360*dihedral/(2.0*3.14159265359)
1868  if(dihedral < 0)dihedral = 360 + dihedral
1869 
1870  end subroutine prg_get_dihedral
1871 
1882  subroutine prg_get_covgraph(sy,nnStruct,nrnnstruct,bml_type,factor,gcov_bml,mdimin,verbose)
1883  implicit none
1884  character(20), intent(in) :: bml_type
1885  integer :: i, j, jj, mdim
1886  integer, intent(in) :: mdimin
1887  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:)
1888  integer, optional, intent(in) :: verbose
1889  real(dp) :: d, dvdw, factor, ra(3),rb(3),rab(3)
1890  real(dp) :: lx, ly, lz
1891  type(bml_matrix_t), intent(inout) :: gcov_bml
1892  type(system_type), intent(in) :: sy
1893  logical(1), allocatable :: ispresent(:)
1894 
1895  if(bml_get_n(gcov_bml).gt.0) call bml_deallocate(gcov_bml)
1896 
1897  if(mdimin > 0)then
1898  mdim = mdimin
1899  else
1900  mdim = sy%nats
1901  endif
1902 
1903  call bml_zero_matrix(bml_type,bml_element_real,kind(1.0),sy%nats,mdim,gcov_bml)
1904  !call bml_zero_matrix(bml_type,bml_element_real,dp,sy%nats,mdim,gcov_bml)
1905 
1906  lx = sy%lattice_vector(1,1)
1907  ly = sy%lattice_vector(2,2)
1908  lz = sy%lattice_vector(3,3)
1909 
1910  if(present(verbose))then
1911  if(verbose >= 1)then
1912  write(*,*)" "
1913  write(*,*)"Building covalency graph ..."
1914  write(*,*)" "
1915  endif
1916  endif
1917 
1918  allocate(ispresent(sy%nats))
1919 
1920  do i=1,sy%nats
1921  ra = sy%coordinate(:,i)
1922  ispresent = .false.
1923  do j=1,nrnnstruct(i)
1924 
1925  if(nrnnstruct(i) <= 1)write(*,*)"WARNING: Atom" ,i,"is desconnected"
1926 
1927  jj = nnstruct(j,i)
1928  if(.not.ispresent(jj))then
1929  rb = sy%coordinate(:,jj)
1930 
1931  rab(1) = modulo((rb(1) - ra(1) + lx/2.0_dp),lx) - lx/2.0_dp
1932  rab(2) = modulo((rb(2) - ra(2) + ly/2.0_dp),ly) - ly/2.0_dp
1933  rab(3) = modulo((rb(3) - ra(3) + lz/2.0_dp),lz) - lz/2.0_dp
1934 
1935  d = prg_norm2(rab)
1936 
1937  if(d == 0.0_dp .and. i .ne. jj)write(*,*)"WARNING: Atom" ,i,"and atom",jj,&
1938  "are on top of each other"
1939  dvdw = factor*(element_covr(sy%atomic_number(i)) + element_covr(sy%atomic_number(jj)))
1940  if(d < dvdw .and. d > 0.0_dp)then
1941  ispresent(jj) = .true.
1942  ! call bml_set_element_new(gcov_bml,i,jj,1.0_dp)
1943  ! call bml_set_element_new(gcov_bml,jj,i,1.0_dp)
1944  call bml_set_element_new(gcov_bml,i,jj,1.0)
1945  ! call bml_set_element_new(gcov_bml,jj,i,1.0)
1946  endif
1947  endif
1948  enddo
1949 
1950  ! call bml_set_element_new(gcov_bml,i,i,1.0_dp)
1951  enddo
1952 
1953  deallocate(ispresent)
1954 
1955  end subroutine prg_get_covgraph
1956 
1957 
1958  subroutine prg_get_covgraph_int(sy,nnStructMindist,nnStruct,nrnnstruct,bml_type,factor,gcov_bml,mdimin,verbose)
1959  implicit none
1960  character(20), intent(in) :: bml_type
1961  integer :: i, j, jj, mdim
1962  integer, intent(in) :: mdimin
1963  integer, intent(in) :: nnStruct(:,:), nrnnstruct(:)
1964  integer, optional, intent(in) :: verbose
1965  real(dp) :: d, dvdw, factor
1966  real(dp), intent(in) :: nnStructMindist(:,:)
1967  type(bml_matrix_t), intent(inout) :: gcov_bml
1968  type(system_type), intent(in) :: sy
1969 
1970  if(bml_get_n(gcov_bml).gt.0) call bml_deallocate(gcov_bml)
1971 
1972  if(mdimin > 0)then
1973  mdim = mdimin
1974  else
1975  mdim = sy%nats
1976  endif
1977 
1978  call bml_zero_matrix(bml_type,bml_element_real,dp,sy%nats,mdim,gcov_bml)
1979 
1980  if(present(verbose))then
1981  if(verbose >= 1)then
1982  write(*,*)" "
1983  write(*,*)"Building covalency graph ..."
1984  write(*,*)" "
1985  endif
1986  endif
1987 
1988  do i=1,sy%nats
1989  do j=1,nrnnstruct(i)
1990  if(nrnnstruct(i) <= 1)write(*,*)"WARNING: Atom" ,i,"is desconnected"
1991  jj = nnstruct(j,i)
1992  d = nnstructmindist(j,i)
1993  if(d == 0.0_dp .and. i .ne. jj)write(*,*)"WARNING: Atom" ,i,"and atom",j,&
1994  "are on top of each other"
1995  dvdw = factor*(element_covr(sy%atomic_number(i)) + element_covr(sy%atomic_number(jj)))
1996  if(d < dvdw .and. d > 0.0_dp)then
1997  call bml_set_element_new(gcov_bml,i,jj,1.0_dp)
1998  call bml_set_element_new(gcov_bml,jj,i,1.0_dp)
1999  endif
2000  enddo
2001  call bml_set_element_new(gcov_bml,i,i,1.0_dp)
2002  enddo
2003 
2004  end subroutine prg_get_covgraph_int
2005 
2006 
2017  subroutine prg_get_covgraph_h(sy,nnStruct,nrnnstruct,rcut,&
2018  graph_h,mdimin,verbose)
2019  implicit none
2020  integer :: i, j, jj, ncount, mdim
2021  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:), mdimin
2022  integer, optional, intent(in) :: verbose
2023  real(dp) :: d, dvdw, ra(3),rb(3),rab(3)
2024  real(dp) :: lx, ly, lz
2025  real(dp), intent(in) :: rcut
2026  integer, allocatable, intent(inout) :: graph_h(:,:)
2027  type(system_type), intent(in) :: sy
2028 
2029 
2030  if(mdimin > 0)then
2031  mdim = mdimin
2032  else
2033  mdim = sy%nats
2034  endif
2035 
2036  if(.not.allocated(graph_h)) allocate(graph_h(mdim,sy%nats))
2037 
2038  if(present(verbose))then
2039  if(verbose >= 1)then
2040  write(*,*)" "
2041  write(*,*)"Building H covalency graph ..."
2042  write(*,*)" "
2043  endif
2044  endif
2045 
2046  graph_h = 0
2047 
2048  lx = sy%lattice_vector(1,1)
2049  ly = sy%lattice_vector(2,2)
2050  lz = sy%lattice_vector(3,3)
2051 
2052  !$omp parallel do default(none) private(i) &
2053  !$omp private(ncount,j,jj,ra,rb,rab,d) &
2054  !$omp shared(sy,nrnnstruct,nnStruct,Lx,Ly,Lz,graph_h,rcut)
2055  do i=1,sy%nats
2056  ncount = 0
2057  do j=1,nrnnstruct(i)
2058  jj = nnstruct(j,i)
2059  ra=sy%coordinate(:,i)
2060  rb=sy%coordinate(:,jj)
2061  rab(1) = modulo((rb(1) - ra(1) + lx/2.0_dp),lx) - lx/2.0_dp
2062  rab(2) = modulo((rb(2) - ra(2) + ly/2.0_dp),ly) - ly/2.0_dp
2063  rab(3) = modulo((rb(3) - ra(3) + lz/2.0_dp),lz) - lz/2.0_dp
2064 
2065  d = prg_norm2(rab)
2066 
2067  if(d.lt.rcut.and.d.gt.0.0_dp)then
2068  ncount=ncount+1
2069  graph_h(ncount,i) = jj
2070  endif
2071  enddo
2072  enddo
2073  !$omp end parallel do
2074 
2075  end subroutine prg_get_covgraph_h
2076 
2087  subroutine prg_get_subsystem(sy,lsize,indices,sbsy,verbose)
2088  implicit none
2089  integer :: i, nsptmp, prev
2090  integer, intent(in) :: indices(:), lsize
2091  integer, optional, intent(in) :: verbose
2092  type(system_type), intent(in) :: sy
2093  type(system_type), intent(inout) :: sbsy
2094 
2095  if(present(verbose))then
2096  if(verbose >= 1)then
2097  write(*,*)" "
2098  write(*,*)"Extracting subsystem ..."
2099  write(*,*)" "
2100  endif
2101  endif
2102 
2103  sbsy%nats = lsize
2104 
2105  if(allocated(sbsy%symbol))then
2106  deallocate(sbsy%symbol)
2107  deallocate(sbsy%coordinate)
2108  deallocate(sbsy%atomic_number)
2109  deallocate(sbsy%velocity)
2110  deallocate(sbsy%force)
2111  deallocate(sbsy%net_charge)
2112  deallocate(sbsy%mass)
2113  deallocate(sbsy%lattice_vector)
2114  deallocate(sbsy%spindex)
2115  endif
2116 
2117  allocate(sbsy%symbol(sbsy%nats))
2118  allocate(sbsy%coordinate(3,sbsy%nats))
2119  allocate(sbsy%atomic_number(sbsy%nats))
2120  allocate(sbsy%velocity(3,sbsy%nats))
2121  allocate(sbsy%force(3,sbsy%nats))
2122  allocate(sbsy%net_charge(sbsy%nats))
2123  allocate(sbsy%mass(sbsy%nats))
2124  allocate(sbsy%lattice_vector(3,3))
2125  allocate(sbsy%spindex(sbsy%nats))
2126 
2127  ! allocate(sbsy%recip_vector(3,3))
2128 
2129  sbsy%lattice_vector = sy%lattice_vector
2130  sbsy%volr = sy%volr
2131  sbsy%volk = sy%volk
2132 
2133  nsptmp = 0
2134  do i=1,sbsy%nats
2135  sbsy%symbol(i) = sy%symbol(indices(i)+1) !Indices from the graph partition start from 0
2136  sbsy%coordinate(:,i) = sy%coordinate(:,indices(i)+1)
2137  sbsy%atomic_number(i) = sy%atomic_number(indices(i)+1)
2138  sbsy%mass(i) = sy%mass(indices(i)+1)
2139  sbsy%spindex(i) = sy%spindex(indices(i)+1)
2140  if(sbsy%spindex(i).gt.nsptmp)then
2141  nsptmp = nsptmp + 1
2142  endif
2143  enddo
2144 
2145  if(nsptmp.lt.sy%nsp)then
2146  write(*,*)"WARNING: nsp = ",nsptmp
2147  write(*,*)"The subsystem contains less species that the system ..."
2148  write(*,*)"A generalization to parts where subsystem contains"
2149  write(*,*)"less species that the system needs to be added ..."
2150  write(*,*)"See prg_get_subsystem in prg_system_mod"
2151  ! stop
2152  endif
2153 
2154  sbsy%nsp = sy%nsp
2155 
2156  if(allocated(sbsy%spatnum))then
2157  deallocate(sbsy%spatnum)
2158  deallocate(sbsy%spmass)
2159  endif
2160 
2161  allocate(sbsy%spatnum(sbsy%nsp))
2162  allocate(sbsy%spmass(sbsy%nsp))
2163 
2164  sbsy%spatnum = sy%spatnum
2165  sbsy%spmass = sy%spmass
2166 
2167  end subroutine prg_get_subsystem
2168 
2169 
2176  subroutine prg_destroy_subsystems(sbsy,verbose)
2177  implicit none
2178  type(system_type), intent(inout) :: sbsy
2179  integer, optional, intent(in) :: verbose
2180  integer :: nparts
2181 
2182  if(present(verbose))then
2183  if(verbose >= 1)then
2184  write(*,*)" "
2185  write(*,*)"Deallocating the multiple subsystem ..."
2186  write(*,*)" "
2187  endif
2188  endif
2189 
2190 
2191  if(allocated(sbsy%symbol))then
2192  deallocate(sbsy%symbol)
2193  deallocate(sbsy%coordinate)
2194  deallocate(sbsy%atomic_number)
2195  deallocate(sbsy%velocity)
2196  deallocate(sbsy%force)
2197  deallocate(sbsy%net_charge)
2198  deallocate(sbsy%mass)
2199  deallocate(sbsy%spindex)
2200  endif
2201 
2202  if(allocated(sbsy%lattice_vector))deallocate(sbsy%lattice_vector)
2203  if(allocated(sbsy%resindex))deallocate(sbsy%resindex)
2204 
2205  if(allocated(sbsy%spatnum))then
2206  deallocate(sbsy%spatnum)
2207  deallocate(sbsy%spmass)
2208  endif
2209 
2210  if(bml_get_n(sbsy%estr%ham) > 0) call bml_deallocate(sbsy%estr%ham)
2211  if(bml_get_n(sbsy%estr%ham0) > 0) call bml_deallocate(sbsy%estr%ham0)
2212  if(bml_get_n(sbsy%estr%oham) > 0) call bml_deallocate(sbsy%estr%oham)
2213  if(bml_get_n(sbsy%estr%over) > 0) call bml_deallocate(sbsy%estr%over)
2214  if(bml_get_n(sbsy%estr%rho) > 0) call bml_deallocate(sbsy%estr%rho)
2215  if(bml_get_n(sbsy%estr%orho) > 0) call bml_deallocate(sbsy%estr%orho)
2216  if(bml_get_n(sbsy%estr%zmat) > 0) call bml_deallocate(sbsy%estr%zmat)
2217 
2218 
2219  if(allocated(sbsy%estr%coul_pot_r)) deallocate(sbsy%estr%coul_pot_r)
2220  if(allocated(sbsy%estr%coul_pot_k)) deallocate(sbsy%estr%coul_pot_k)
2221  if(allocated(sbsy%estr%skforce)) deallocate(sbsy%estr%skforce)
2222  if(allocated(sbsy%estr%fpul)) deallocate(sbsy%estr%fpul)
2223  if(allocated(sbsy%estr%fscoul)) deallocate(sbsy%estr%fscoul)
2224  if(allocated(sbsy%estr%hindex)) deallocate(sbsy%estr%hindex)
2225 
2226 
2227  end subroutine prg_destroy_subsystems
2228 
2240  subroutine prg_molpartition(sy,npart,nnStructMindist,nnStruct,nrnnstruct,hetatm,gp,verbose)
2241  implicit none
2242  character(2), intent(in) :: hetatm
2243  integer :: cnt, i, j, jj
2244  integer :: maxparts, nmax
2245  integer, allocatable :: core_count(:), cores(:,:), partof(:)
2246  integer, intent(in) :: nnstruct(:,:), nrnnstruct(:)
2247  integer, intent(inout) :: npart
2248  integer, optional, intent(inout) :: verbose
2249  logical, allocatable :: inpart(:)
2250  real(dp) :: d, dvdw
2251  real(dp), intent(in) :: nnstructmindist(:,:)
2252  type(graph_partitioning_t), intent(inout) :: gp
2253  type(system_type), intent(in) :: sy
2254 
2255  if(present(verbose))then
2256  if(verbose >= 1)then
2257  write(*,*)" "; write(*,*)"Partitioning by molecule ..."; write(*,*)" "
2258  endif
2259  endif
2260 
2261  if(allocated(gp%nnodesInPart))then
2262  call prg_destroygraphpartitioning(gp)
2263  endif
2264 
2265  nmax = sy%nats
2266  maxparts = sy%nats
2267 
2268  allocate(core_count(maxparts))
2269 
2270  allocate(cores(nmax,maxparts))
2271  allocate(inpart(sy%nats))
2272  allocate(partof(sy%nats))
2273 
2274  cores = -1
2275  inpart = .false.
2276  partof=-1
2277  core_count = 0
2278  npart = 0
2279  do i=1,sy%nats
2280  if(.not.inpart(i).and.trim(sy%symbol(i)) == trim(hetatm))then
2281  npart = npart + 1
2282  core_count(npart) = core_count(npart) + 1
2283  cores(core_count(npart),npart) = i
2284  partof(i) = npart
2285  do j=1,nrnnstruct(i)
2286  jj = nnstruct(j,i)
2287  d = nnstructmindist(j,i)
2288  dvdw = 1.3_dp
2289  if(d.lt.dvdw.and.d.gt.0.001_dp)then
2290  if(.not.inpart(jj))then
2291  core_count(npart) = core_count(npart) + 1
2292  cores(core_count(npart),npart) = jj
2293  inpart(jj) = .true.
2294  partof(jj) = npart
2295  endif
2296  endif
2297  enddo
2298  endif
2299  enddo
2300 
2301  call prg_initgraphpartitioning(gp, "molecules", npart, sy%nats, sy% nats)
2302 
2303  ! Initialize and fill up subgraph structure
2304  ! Assign node ids (mapped to orbitals as rows) to each node in each
2305  do i = 1, npart
2306  gp%nnodesInPartAll(i) = core_count(i)
2307  call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
2308  allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
2309  gp%nnodesInPart(i) = core_count(i)
2310  enddo
2311 
2312  !Assign node ids to sgraph
2313  do i=1, npart
2314  do j=1,core_count(i)
2315  gp%sgraph(i)%nodeInPart(j) = cores(j,i) - 1
2316  end do
2317  enddo
2318 
2319  end subroutine prg_molpartition
2320 
2329  subroutine prg_get_partial_atomgraph(rho_bml,hindex,gch_bml,threshold,verbose)
2330  implicit none
2331  character(20) :: bml_type
2332  integer :: i, ii, j, jj
2333  integer :: nch, norbs
2334  integer, intent(in) :: hindex(:,:)
2335  integer, optional, intent(in) :: verbose
2336  logical :: connection
2337  logical, allocatable :: iconnectedtoj(:)
2338  real(dp), allocatable :: row(:), rowat(:)
2339  real(dp), intent(in) :: threshold
2340  type(bml_matrix_t), intent(in) :: rho_bml
2341  type(bml_matrix_t), intent(inout) :: gch_bml
2342 
2343  norbs = bml_get_n(rho_bml)
2344  bml_type = bml_get_type(rho_bml)
2345  nch = size(hindex,dim=2)
2346  allocate(rowat(nch))
2347  allocate(row(norbs))
2348  allocate(iconnectedtoj(nch))
2349 
2350  if(bml_get_n(gch_bml) > 0) call bml_deallocate(gch_bml)
2351  call bml_zero_matrix(bml_type,bml_element_real,dp,nch,nch,gch_bml)
2352 
2353  do i = 1, nch
2354  rowat = 0.0_dp
2355  do ii = hindex(1,i),hindex(2,i) ! i,j block
2356  call bml_get_row(rho_bml,ii,row)
2357  iconnectedtoj = .false.
2358  do j = 1, nch
2359  connection = .false.
2360 
2361  if(.not.iconnectedtoj(j).and.rowat(j) == 0)then
2362  do jj = hindex(1,j),hindex(2,j)
2363  if(abs(row(jj)) > threshold)then
2364  connection = .true. !We exit if there is a connection.
2365  exit
2366  endif
2367  enddo
2368  if(connection) iconnectedtoj(j) = .true.
2369  endif
2370  if(connection) rowat(j) = 1.0_dp
2371  enddo
2372 
2373  enddo
2374  call bml_set_row(gch_bml,i,rowat)
2375  enddo
2376 
2377  deallocate(rowat)
2378  deallocate(row)
2379  deallocate(iconnectedtoj)
2380 
2381  end subroutine prg_get_partial_atomgraph
2382 
2383 
2395  subroutine prg_collect_graph_p(rho_bml,nc,nats,hindex,chindex,graph_p,threshold,mdimin,verbose)
2396  implicit none
2397  character(20) :: bml_type
2398  integer :: i, ifull, ii, j
2399  integer :: jfull, jj, nch, ncounti
2400  integer :: norbs, mdim
2401  logical(1), allocatable :: rowatfull(:)
2402  integer, allocatable, intent(inout) :: graph_p(:,:)
2403  integer, intent(in) :: chindex(:), hindex(:,:), nats, nc
2404  integer, intent(in) :: mdimin
2405  integer, optional, intent(in) :: verbose
2406  logical :: connection
2407  logical, allocatable :: iconnectedtoj(:)
2408  real(dp), allocatable :: row(:)
2409  real(dp), intent(in) :: threshold
2410  type(bml_matrix_t), intent(in) :: rho_bml
2411 
2412  norbs = bml_get_n(rho_bml)
2413  bml_type = bml_get_type(rho_bml)
2414  nch = size(hindex,dim=2)
2415  allocate(rowatfull(nats))
2416  allocate(row(norbs))
2417  allocate(iconnectedtoj(nch))
2418 
2419  if(mdimin > 0)then
2420  mdim = mdimin
2421  else
2422  mdim = nats
2423  endif
2424 
2425  if(.not.allocated(graph_p))then
2426  allocate(graph_p(mdim,nats))
2427  graph_p = 0
2428  endif
2429 
2430  !$omp parallel do default(none) private(i) &
2431  !$omp private(ncounti,rowatfull,ii,j,jfull,ifull,iconnectedtoj) &
2432  !$omp private(row,connection) &
2433  !$omp shared(graph_p,nc,chindex,hindex,rho_bml,nch,threshold)
2434  do i = 1, nc
2435 
2436  ifull = chindex(i) + 1 !Map it to the full system
2437  rowatfull = .false.
2438  ii=1
2439  ncounti = 0
2440  do while (graph_p(ii,ifull).ne.0) !Unfolding the connections of ifull
2441  rowatfull(graph_p(ii,ifull)) = .true.
2442  ncounti = ncounti + 1
2443  ii = ii+1
2444  enddo
2445 
2446  do ii = hindex(1,i),hindex(2,i) ! i,j block
2447  call bml_get_row(rho_bml,ii,row)
2448  iconnectedtoj = .false.
2449  !
2450  do j = 1, nch
2451  jfull = chindex(j) + 1 !Cause the count starts form 0
2452  connection = .false.
2453  if(.not.iconnectedtoj(j).and..not.rowatfull(jfull))then
2454  do jj = hindex(1,j),hindex(2,j)
2455  if(abs(row(jj)) > threshold)then
2456  connection = .true. !We exit if there is a connection.
2457  exit
2458  endif
2459  enddo
2460  if(connection) iconnectedtoj(j) = .true.
2461  endif
2462  if(connection)then !I conected to j >>> ifull connected to jfull
2463  !Map back to the full system graph.
2464  rowatfull(jfull) = .true.
2465  ncounti = ncounti + 1
2466  graph_p(ncounti,ifull) = jfull
2467  endif
2468  enddo
2469  !
2470  enddo
2471 
2472  enddo
2473  !$omp end parallel do
2474 
2475  deallocate(rowatfull)
2476  deallocate(iconnectedtoj)
2477  deallocate(row)
2478 
2479  end subroutine prg_collect_graph_p
2480 
2481 
2487  subroutine prg_merge_graph(graph_p,graph_h)
2488  implicit none
2489  integer :: i, ii, j, nats
2490  integer :: ncounti, maxnz
2491  integer, intent(inout) :: graph_h(:,:)
2492  integer, intent(inout) :: graph_p(:,:)
2493  logical(1), allocatable :: rowpatfull(:)
2494 
2495  nats = size(graph_p,dim=2)
2496  maxnz = size(graph_p,dim=1)
2497 
2498  allocate(rowpatfull(nats))
2499 
2500  !$omp parallel do default(none) private(i) &
2501  !$omp private(ncounti,rowpatfull,ii,j) &
2502  !$omp shared(graph_p,graph_h,nats,maxnz)
2503  do i = 1, nats
2504  ncounti = 0
2505  rowpatfull = .false.
2506 
2507  do ii = 1,maxnz
2508  if(graph_p(ii,i) == 0) exit
2509  rowpatfull(graph_p(ii,i)) = .true.
2510  ncounti = ncounti + 1
2511  enddo
2512 
2513  do ii = 1,maxnz
2514  j = graph_h(ii,i)
2515  if(j==0)exit
2516  if(.not.rowpatfull(j))then
2517  ncounti = ncounti + 1
2518  graph_p(ncounti,i) = j
2519  endif
2520  enddo
2521 
2522  enddo
2523  !$omp end parallel do
2524 
2525  deallocate(rowpatfull)
2526 
2527 
2528  end subroutine prg_merge_graph
2529 
2530 
2538  subroutine prg_merge_graph_adj(graph_p,graph_h,xadj,adjncy)
2539  implicit none
2540  integer :: i, ii, j, nats
2541  integer :: ncounti, ncountot
2542  integer, allocatable, intent(inout) :: adjncy(:), graph_h(:,:), graph_p(:,:), xadj(:)
2543  logical(1), allocatable :: rowpatfull(:)
2544 
2545  nats = size(graph_p,dim=2)
2546  allocate(rowpatfull(nats))
2547  allocate(xadj(nats+1))
2548  allocate(adjncy(nats*nats))
2549 
2550  xadj = 0.0
2551  adjncy = 0.0
2552 
2553  !$omp parallel do default(none) private(i) &
2554  !$omp private(ncounti,rowpatfull,ii,j,ncountot) &
2555  !$omp shared(graph_p,graph_h,nats)
2556  do i = 1, nats
2557  ncounti = 0
2558  rowpatfull = .false.
2559  ii = 1
2560  do while (graph_p(ii,i).ne.0) !Unfolding the connections of i
2561  rowpatfull(graph_p(ii,i)) = .true.
2562  ncounti = ncounti + 1
2563  ii = ii+1
2564  enddo
2565 
2566  ii = 1
2567  do while (graph_h(ii,i).ne.0) !Unfolding the connections of i
2568  j = graph_h(ii,i)
2569  if(.not.rowpatfull(j))then
2570  ncounti = ncounti + 1
2571  ncountot = ncountot + 1
2572  graph_p(ncounti,i) = j
2573  endif
2574  ii = ii+1
2575  enddo
2576  enddo
2577  !$omp end parallel do
2578 
2579 
2580  ncountot = 0
2581 
2582  !$omp parallel do default(none) private(i) &
2583  !$omp private(ncounti,ii,ncountot) &
2584  !$omp shared(graph_p,nats,xadj,adjncy)
2585  do i = 1, nats
2586  ncounti = 0
2587  ii = 1
2588  do while (graph_p(ii,i).gt.0)
2589  ncounti = ncounti + 1
2590  ncountot = ncountot + 1
2591  adjncy(ncountot) = graph_p(ii,i)
2592  ii = ii+1
2593  enddo
2594  xadj(i) = ncountot - ncounti + 1
2595  enddo
2596  !$omp end parallel do
2597 
2598  xadj(nats+1) = ncountot + 1
2599 
2600  deallocate(rowpatfull)
2601  deallocate(graph_p)
2602 
2603  end subroutine prg_merge_graph_adj
2604 
2612  subroutine prg_adj2bml(xadj,adjncy,bml_type,g_bml)
2613  implicit none
2614  character(20), intent(in) :: bml_type
2615  integer :: i, ii, j, nats
2616  integer, intent(in) :: adjncy(:), xadj(:)
2617  real(dp), allocatable :: row(:)
2618  type(bml_matrix_t), intent(inout) :: g_bml
2619 
2620  nats = size(xadj,dim=1) - 1
2621  if(bml_get_n(g_bml).gt.0) call bml_deallocate(g_bml)
2622  call bml_zero_matrix(bml_type,bml_element_real,dp,nats,nats,g_bml)
2623  allocate(row(nats))
2624 
2625  !$omp parallel do default(none) private(i) &
2626  !$omp private(j,row) &
2627  !$omp shared(nats,g_bml,xadj,adjncy)
2628  do i = 1, nats
2629  row = 0.0_dp
2630  do j = xadj(i), xadj(i+1) - 1
2631  row(adjncy(j)) = 1.0_dp
2632  enddo
2633  call bml_set_row(g_bml,i,row)
2634  enddo
2635 
2636  deallocate(row)
2637 
2638  end subroutine prg_adj2bml
2639 
2646  subroutine prg_graph2bml(graph,bml_type,g_bml)
2647  implicit none
2648  character(20), intent(in) :: bml_type
2649  integer :: i, ii, nats, mdim
2650  integer, allocatable, intent(inout) :: graph(:,:)
2651  !real(dp), allocatable :: row(:)
2652  real(kind(1.0)), allocatable :: row(:)
2653  type(bml_matrix_t), intent(inout) :: g_bml
2654 
2655  nats = size(graph,dim=2)
2656  ! if(bml_get_N(g_bml).GT.0) call bml_deallocate(g_bml)
2657  ! call bml_zero_matrix(bml_type,bml_element_real,kind(1.0),nats,nats,g_bml)
2658  allocate(row(nats))
2659 
2660  mdim = bml_get_m(g_bml)
2661 
2662  !$omp parallel do default(none) private(i) &
2663  !$omp private(ii,row) &
2664  !$omp shared(nats,g_bml,graph,mdim)
2665  do i = 1, nats
2666  ii = 1
2667  row = 0.0
2668  do while (graph(ii,i).gt.0)
2669  row(graph(ii,i)) = 1.0
2670  ! call bml_set_element_new(g_bml,i,graph(ii,i),1.0)
2671  ! call bml_set_element_new(g_bml,graph(ii,i),i,1.0)
2672  ii = ii+1
2673  enddo
2674  call bml_set_row(g_bml,i,row,0.5_dp)
2675  enddo
2676  !$omp end parallel do
2677 
2678  deallocate(graph)
2679  deallocate(row)
2680 
2681  end subroutine prg_graph2bml
2682 
2683 
2689  subroutine prg_graph2vector(graph,vector,maxnz)
2690  implicit none
2691  integer, intent(inout) :: graph(:,:)
2692  integer, allocatable :: vector(:)
2693  integer :: i, j, maxnz, nats, ncount
2694 
2695  nats = size(graph,dim=2)
2696  allocate(vector(nats*maxnz))
2697 
2698  ncount = 0
2699 
2700  !$omp parallel do default(none) private(i) &
2701  !$omp private(ncount,j) &
2702  !$omp shared(nats,maxnz,vector,graph)
2703  do i = 1, nats
2704  do j = 1, maxnz
2705  ncount = (i-1)*maxnz + j !ncount + 1
2706  vector(ncount) = graph(j,i)
2707  enddo
2708  enddo
2709  !$omp end parallel do
2710 
2711  end subroutine prg_graph2vector
2712 
2718  subroutine prg_vector2graph(vector,graph,maxnz)
2719  implicit none
2720  integer, intent(inout) :: graph(:,:)
2721  integer, allocatable, intent(inout) :: vector(:)
2722  integer :: i, j, maxnz, nats, ncount
2723 
2724  nats = size(graph,dim=2)
2725  ncount = 0
2726 
2727  !$omp parallel do default(none) private(i) &
2728  !$omp private(ncount,j) &
2729  !$omp shared(nats,maxnz,vector,graph)
2730  do i = 1, nats
2731  do j = 1, maxnz
2732  ncount = (i-1)*maxnz + j !ncount + 1
2733  graph(j,i) = vector(ncount)
2734  enddo
2735  enddo
2736  !$omp end parallel do
2737 
2738  deallocate(vector)
2739 
2740  end subroutine prg_vector2graph
2741 
2746  subroutine prg_sortadj(xadj, adjncy)
2747  implicit none
2748  integer, intent(inout) :: xadj(:)
2749  integer, allocatable, intent(inout) :: adjncy(:)
2750  integer :: i, j, n, l, first, last, nx, irank, mintmp
2751  integer, allocatable :: tmpvect(:), newadjncy(:)
2752 
2753  nx = size(xadj,dim=1)
2754  allocate(tmpvect(nx))
2755  allocate(newadjncy(xadj(nx) - 1))
2756 
2757  !$omp parallel do default(none) private(i) &
2758  !$omp private(first,last,j,l,tmpvect,N,mintmp) &
2759  !$omp shared(nx,adjncy,xadj,newadjncy)
2760  do i = 1, nx-1
2761  first = xadj(i)
2762  last = xadj(i+1) - 1
2763  write(*,*)first,last
2764  n = last - first + 1
2765  tmpvect(1:n) = 0
2766  tmpvect(1:n) = adjncy(first:last)
2767 
2768  do j = 1,n
2769  do l = j+1,n
2770  if(tmpvect(j) >= tmpvect(l).and.tmpvect(l).ne.0)then
2771  mintmp = tmpvect(l)
2772  tmpvect(l) = tmpvect(j)
2773  tmpvect(j) = mintmp
2774  endif
2775  enddo
2776  enddo
2777 
2778  adjncy(first:last) = tmpvect(1:n)
2779  newadjncy(first:last) = tmpvect(1:n)
2780 
2781  enddo
2782  !$omp end parallel do
2783  deallocate(tmpvect)
2784  deallocate(adjncy)
2785 
2786  allocate(adjncy(xadj(nx) - 1)) !Reallocation for adjusting size
2787 
2788  adjncy = newadjncy
2789 
2790  deallocate(newadjncy)
2791 
2792  end subroutine prg_sortadj
2793 
2794 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:175
prg_graph_mod::prg_initsubgraph
subroutine, public prg_initsubgraph(sg, pnum, hsize)
Initialize subgraph.
Definition: prg_graph_mod.F90:143
prg_graph_mod::prg_destroygraphpartitioning
subroutine, public prg_destroygraphpartitioning(gp)
Destroy graph partitioning.
Definition: prg_graph_mod.F90:263
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:2088
prg_system_mod::prg_adj2bml
subroutine, public prg_adj2bml(xadj, adjncy, bml_type, g_bml)
prg_adj2bml
Definition: prg_system_mod.F90:2613
prg_system_mod::prg_destroy_subsystems
subroutine, public prg_destroy_subsystems(sbsy, verbose)
Destroy allocated subsystem.
Definition: prg_system_mod.F90:2177
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:1229
prg_system_mod::prg_translatetogeomcandfoldtobox
subroutine, public prg_translatetogeomcandfoldtobox(coords, lattice_vectors, origin)
Translate to geometric center.
Definition: prg_system_mod.F90:1505
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:1601
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:1828
prg_system_mod::prg_graph2bml
subroutine, public prg_graph2bml(graph, bml_type, g_bml)
Graph2bml.
Definition: prg_system_mod.F90:2647
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:950
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:1959
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:2241
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:1183
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:1698
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:2539
prg_system_mod::prg_graph2vector
subroutine, public prg_graph2vector(graph, vector, maxnz)
Vectorize graph.
Definition: prg_system_mod.F90:2690
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:1883
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:1073
prg_system_mod::prg_vector2graph
subroutine, public prg_vector2graph(vector, graph, maxnz)
Back to graph.
Definition: prg_system_mod.F90:2719
prg_system_mod::prg_translateandfoldtobox
subroutine, public prg_translateandfoldtobox(coords, lattice_vectors, origin, verbose)
Translate and fold to box.
Definition: prg_system_mod.F90:1365
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:1460
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:705
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:2019
prg_system_mod::prg_get_distancematrix
subroutine, public prg_get_distancematrix(coords, dmat)
Get the distance matrix.
Definition: prg_system_mod.F90:1340
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:245
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:2330
prg_system_mod::prg_destroy_system
subroutine, public prg_destroy_system(sy)
Deallocates all the arrays within a system.
Definition: prg_system_mod.F90:678
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:268
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:1418
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:1271
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:2488
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:1544
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:738
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:1779
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:2747
prg_system_mod::prg_get_origin
subroutine, public prg_get_origin(coords, origin)
Get the origin of the coordinates.
Definition: prg_system_mod.F90:1307
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:2396
prg_graph_mod::graph_partitioning_t
Trace per iteration.
Definition: prg_graph_mod.F90:57
prg_system_mod::estruct_type
Electronic structure type.
Definition: prg_system_mod.F90:20
prg_system_mod::system_type
System type.
Definition: prg_system_mod.F90:107