LCOV - code coverage report
Current view: top level - io - relax_io.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 35 96 36.5 %
Date: 2019-09-08 04:53:50 Functions: 4 5 80.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_relaxio
       8             :   !This module handles IO to the relax.xml file
       9             :   !The writing is done directly to relax.xml
      10             :   !The reading uses the libxml interface to inp.xml. Hence the relax.xml has to be included here.
      11             :   USE m_judft
      12             :   IMPLICIT NONE
      13             :   PRIVATE
      14             :   PUBLIC :: read_relax,write_relax,apply_displacements,read_displacements
      15             : CONTAINS
      16           1 :   SUBROUTINE write_relax(positions,forces,energies,displace)
      17             :     REAL,INTENT(in):: positions(:,:,:)
      18             :     REAL,INTENT(in):: forces(:,:,:)
      19             :     REAL,INTENT(in):: energies(:)
      20             :     REAL,INTENT(in):: displace(:,:)
      21             : 
      22             :     INTEGER :: no_steps,n,ntype,step
      23           1 :     No_steps=SIZE(positions,3)
      24           1 :     ntype=SIZE(positions,2)
      25             :     IF (ntype.NE.SIZE(forces,2).OR.ntype.NE.SIZE(displace,2).OR.&
      26           1 :          no_steps.NE.SIZE(forces,3).OR.no_steps.NE.SIZE(energies))THEN
      27           0 :        CALL judft_error("BUG in relax_io")
      28             :     ENDIF
      29           1 :     OPEN(765,file="relax.xml",status="replace")
      30           1 :     WRITE(765,*) "<!-- Attention, absolute coordinates used here -->"
      31           1 :     WRITE(765,*) "<relaxation>"
      32             :     !write current set of displacements
      33           1 :     WRITE(765,*) "  <displacements>"
      34           3 :     DO n=1,SIZE(displace,2)
      35             :        WRITE(765,"(a,3(f15.10,1x),a)") &
      36           3 :             '    <displace>',displace(:,n),'</displace>'
      37             :     END DO
      38           1 :     WRITE(765,"(a)") '  </displacements>'
      39             : 
      40             :     !Write all known old positions,forces and energies
      41           1 :     WRITE(765,*) "  <relaxation-history>"
      42           2 :     DO step=1,no_steps
      43           1 :        WRITE(765,"(a,f20.10,a)") '    <step energy="',energies(step),'">'
      44           3 :        DO n=1,ntype
      45             :           WRITE(765,"(a,6(f15.10,1x),a)") &
      46           3 :                '      <posforce>',positions(:,n,step),forces(:,n,step),'</posforce>'
      47             :        END DO
      48           2 :        WRITE(765,"(a)") '    </step>'
      49             :     ENDDO
      50           1 :     WRITE(765,*) "  </relaxation-history>"
      51           1 :     WRITE(765,*) "</relaxation>"
      52           1 :     CLOSE(765)
      53           1 :   END SUBROUTINE write_relax
      54             : 
      55           1 :   SUBROUTINE read_relax(positions,forces,energies)
      56             :     USE m_xmlIntWrapFort 
      57             :     USE m_calculator
      58             :     REAL,INTENT(INOUT),ALLOCATABLE:: positions(:,:,:)
      59             :     REAL,INTENT(INOUT),ALLOCATABLE:: forces(:,:,:)
      60             :     REAL,INTENT(INOUT),ALLOCATABLE:: energies(:)
      61             : 
      62           1 :     REAL,ALLOCATABLE::rtmp(:,:,:)
      63             :     INTEGER:: no_steps
      64             :     INTEGER:: ntype,step,n
      65             :     CHARACTER(len=100):: path,p,str
      66           1 :     no_steps=xmlGetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
      67           1 :     ntype=SIZE(positions,2)
      68           1 :     IF (no_steps==0) THEN
      69           1 :        IF (.NOT.ALLOCATED(positions)) ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
      70             :        RETURN
      71             :     END IF
      72           0 :     IF (ALLOCATED(positions)) THEN 
      73             :        !Assume that we got already a set of positions, forces, energy and extend that list
      74           0 :        rtmp=positions
      75           0 :        DEALLOCATE(positions)
      76           0 :        ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
      77           0 :        positions(:,:,no_steps+1:)=rtmp
      78           0 :        rtmp=forces
      79           0 :        DEALLOCATE(forces)
      80           0 :        ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
      81           0 :        forces(:,:,no_steps+1:)=rtmp
      82           0 :        rtmp(1,1,:)=energies
      83           0 :        DEALLOCATE(energies)
      84           0 :        ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
      85           0 :        energies(no_steps+1:)=rtmp(1,1,:)
      86             :     ELSE
      87           0 :        ALLOCATE(positions(3,ntype,no_steps))
      88           0 :        ALLOCATE(forces(3,ntype,no_steps))
      89           0 :        ALLOCATE(energies(no_steps))
      90             :     END IF
      91           0 :     DO step=1,no_steps
      92           0 :        WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/relaxation-history/step[',step,']'
      93           0 :        energies(step)=evaluateFirstOnly(xmlGetAttributeValue(TRIM(path)//"/@energy"))
      94           0 :        DO n=1,ntype
      95           0 :           WRITE(p,"(a,a,i0,a)") TRIM(path),"/posforce[",n,"]"
      96           0 :           str=xmlGetAttributeValue(p)
      97           0 :           positions(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
      98           0 :           Forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
      99             :        ENDDO
     100             :     END DO
     101           0 :   END SUBROUTINE read_relax
     102             : 
     103             : 
     104          25 :   SUBROUTINE read_displacements(atoms,disp)
     105             :     USE m_xmlIntWrapFort 
     106             :     USE m_calculator
     107             :     USE m_types
     108             :     TYPE(t_atoms),INTENT(in)::atoms
     109             :     REAL,INTENT(out)::disp(:,:)
     110             :     CHARACTER(len=50):: path,str
     111             :     INTEGER :: n
     112          82 :     disp=0.0
     113          50 :     IF (xmlGetNumberOfNodes('/fleurInput/relaxation/displacements')==0) RETURN
     114             :     !read displacements and apply to positions
     115           0 :     IF (atoms%ntype.NE.xmlGetNumberOfNodes('/fleurInput/relaxation/displacements/displace')) CALL judft_error("Wrong number of displacements in relaxation")
     116           0 :     DO n=1,atoms%ntype
     117           0 :        WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/displacements/displace[',n,']'
     118           0 :        str=xmlGetAttributeValue(path)
     119           0 :        disp(:,n)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
     120             :     END DO
     121             :   END SUBROUTINE read_displacements
     122             : 
     123          24 :   SUBROUTINE apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
     124             :     USE m_types
     125             :     USE m_chkmt
     126             :     USE m_constants
     127             :     USE m_mapatom
     128             :     TYPE(t_input),INTENT(IN)   :: input
     129             :     TYPE(t_vacuum),INTENT(IN)  :: vacuum
     130             :     TYPE(t_cell),INTENT(IN)    :: cell
     131             :     TYPE(t_oneD),INTENT(IN)    :: oneD
     132             :     TYPE(t_sym),INTENT(INOUT)  :: sym
     133             :     TYPE(t_noco),INTENT(IN)    :: noco
     134             : 
     135             :     TYPE(t_atoms),INTENT(INOUT):: atoms
     136             : 
     137             : 
     138             :     INTEGER:: n,indx(2)
     139          24 :     REAL   :: disp(3,atoms%ntype),disp_all(3,atoms%nat),taual0(3,atoms%nat),overlap(0:atoms%ntype,atoms%ntype)
     140             : 
     141          24 :     CALL read_displacements(atoms,disp)
     142             :     !change displacements to internal coordinates
     143          24 :     disp=MATMUL(cell%bmat,disp)/tpi_const
     144          79 :     IF (ALL(ABS(disp)<1E-8)) RETURN
     145             :     !Fist make sure original MT spheres do not overlap
     146           0 :     CALL chkmt(atoms,input,vacuum,cell,oneD,.TRUE.,overlap=overlap)
     147           0 :     IF(ANY(overlap>0.0)) CALL judft_error("Overlapping MT-spheres in relaxation before even applying any displacement",hint="You messed up your setup")
     148             : 
     149           0 :     taual0=atoms%taual !Store original positions
     150             : 
     151             :     !Now check for overlapping mt-spheres
     152           0 :     overlap=1.0
     153           0 :     DO WHILE(ANY(overlap>1E-10))
     154           0 :        atoms%taual=taual0  
     155           0 :        CALL rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
     156           0 :        atoms%taual=taual0+disp_all
     157           0 :        atoms%pos=MATMUL(cell%amat,atoms%taual)
     158           0 :        CALL chkmt(atoms,input,vacuum,cell,oneD,.TRUE.,overlap=overlap)
     159           0 :        IF (ANY(overlap>0.0)) THEN
     160           0 :           IF (ANY(overlap(0,:)>1E-10)) CALL judft_error("Atom spills out into vacuum during relaxation")
     161           0 :           indx=MAXLOC(overlap(1:,:)) !the two indices with the most overlap
     162             :           !Try only 90% of displacement
     163           0 :           disp(:,indx(1))=disp(:,indx(1))*0.9 
     164           0 :           disp(:,indx(2))=disp(:,indx(2))*0.9
     165           0 :           WRITE(*,*) "Attention: Overlapping MT-spheres. Reduced displacement by 10%"
     166           0 :           WRITE(*,*) indx,overlap(indx(1),indx(2))
     167             :        END IF
     168             :     END DO
     169             : 
     170           0 :     CALL mapatom(sym,atoms,cell,input,noco)
     171             : 
     172           0 :     WRITE(6,*) "Atomic positions including displacements:"
     173           0 :     DO n=1,atoms%nat
     174           0 :        WRITE(6,"(i4,6(1x,f12.5))") n,atoms%taual(:,n),atoms%pos(:,n)
     175             :     ENDDO
     176             : 
     177             :   END SUBROUTINE apply_displacements
     178             : 
     179           0 :   SUBROUTINE rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
     180             :     USE m_types
     181             :     REAL,INTENT(in)          :: disp(:,:)
     182             :     TYPE(t_atoms),INTENT(in) :: atoms
     183             :     TYPE(t_cell),INTENT(IN)  :: cell
     184             :     TYPE(t_sym),INTENT(IN)   :: sym
     185             :     REAL,INTENT(out)         :: disp_all(:,:)
     186             : 
     187             :     INTEGER:: n,na,jop
     188             :     REAL   :: tau0(3),tau0_rot(3),tau_rot(3)
     189             : 
     190             : 
     191           0 :     DO n=1,atoms%ntype
     192           0 :        tau0=atoms%taual(:,n)
     193           0 :        DO na=SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
     194           0 :           jop = sym%invtab(atoms%ngopr(na))
     195           0 :           tau0_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0)+sym%tau(:,jop) !translation will cancel, included for clarity
     196           0 :           tau_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0+disp(:,n))+sym%tau(:,jop)
     197           0 :           disp_all(:,na)=tau_rot-tau0_rot
     198             :        END DO
     199             :     END DO
     200           0 :   END SUBROUTINE rotate_to_all_sites
     201             : END MODULE m_relaxio

Generated by: LCOV version 1.13