TABLE OF CONTENTS


ABINIT/critics [ Functions ]

[ Top ] [ Functions ]

NAME

 critics

FUNCTION

 Search for critical points starting between
    atom inxat and its neighbors.

COPYRIGHT

 Copyright (C) 2002-2017 ABINIT group (PCasek,FF,XG,MM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  aim_dtset= the structured entity containing all input variables
  dstmax=maximum distance to search for neighbors
  stwo, sthree, sfour: logical switches (TRUE/FALSE) indicating
                          to search CP starting in the middle point
                          of two, three or four atoms. One of these
                          atoms is inxat.

OUTPUT

  (see side effects)

SIDE EFFECTS

  This routines acts primarily on the data contained in the aim_prom module

 WARNING
 This file does not follow the ABINIT coding rules (yet)

PARENTS

      drvaim

CHILDREN

      bschg1,critic,sort_dp,vgh_rho

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 
 49 subroutine  critics(aim_dtset,inxat,stwo,sthree,sfour,dstmax)
 50 
 51  use defs_basis
 52  use defs_aimprom
 53  use defs_parameters
 54  use defs_abitypes
 55  use m_sort
 56  use m_profiling_abi
 57 
 58 !This section has been created automatically by the script Abilint (TD).
 59 !Do not modify the following lines by hand.
 60 #undef ABI_FUNC
 61 #define ABI_FUNC 'critics'
 62  use interfaces_63_bader, except_this_one => critics
 63 !End of the abilint section
 64 
 65  implicit none
 66 
 67 !Arguments ------------------------------------
 68 !scalars
 69  integer,intent(in) :: inxat
 70  real(dp),intent(in) :: dstmax
 71  logical,intent(in) :: sfour,sthree,stwo
 72 !no_abirules
 73  type(aim_dataset_type), intent(in) :: aim_dtset
 74 
 75 !Local variables ------------------------------
 76 !scalars
 77  integer :: i1,i2,i3,iat,ii,ipos,ires,jii,jj,kjj,kk,ll,n1,n2,n3,nb,nc
 78  integer :: nshell
 79  real(dp) :: chg,dif1,dif2,diff,dist,olddist,rr
 80 ! real(dp) :: ss,uu
 81  logical :: found,inter
 82 !arrays
 83  integer :: ibat(nnpos*natom),inat(nnpos*natom),ipibat(nnpos*natom)
 84  integer :: nnat(nnpos*natom),nr(nnpos*natom)
 85  real(dp) :: dif(3),dists(nnpos*natom),ev(3),grho(3),hrho(3,3)
 86  real(dp) :: pom(3),v1(3),v2(3),v3(3),v4(3),vi(3),vt(3),zz(3,3)
 87 
 88 !************************************************************************
 89  vi(:)=xatm(:,inxat)
 90 
 91  nc=0
 92  do jii=1,nnpos
 93    do kjj=1,natom
 94      dist=0._dp
 95      dif(:)=xatm(:,inxat)-xatm(:,kjj)-atp(:,jii)
 96 
 97 !    do ii=1,3
 98 !    dif(ii)=xatm(ii,inxat)-xatm(ii,kjj)-atp(ii,jii)
 99 !    end do
100      dist=vnorm(dif,0)
101      if (.not.((dist>dstmax).or.(dist<0.001))) then
102        nc=nc+1
103        dists(nc)=dist
104        nnat(nc)=kjj
105        inat(nc)=jii
106      end if
107    end do
108  end do
109  do n1=1,nc
110    nr(n1)=n1
111  end do
112  call sort_dp(nc,dists,nr,tol14)
113  nb=0
114  olddist=0._dp
115  nshell=0
116 !write(std_out,*) ':ORIAT ', (xatm(ii,inxat),ii=1,3)
117  do n1=1,nc
118    n2=nr(n1)
119    n3=nnat(n2)
120    if (dists(n1)<(2*dists(1))) then
121      if ((dists(n1)-olddist)>aim_dlimit) then
122        nshell=nshell+1
123        olddist=dists(n1)
124        if (nshell==5) exit
125      end if
126      nb=nb+1
127      ibat(nb)=n3
128      ipibat(nb)=inat(n2)
129      write(std_out,*) ':NEIG ',inxat,n3,inat(n2),dists(n1)
130 !    write(std_out,*) ':POSAT',(xatm(ii,ibat(nb))+atp(ii,ipibat(nb)),ii=1,3)
131    else
132      exit
133    end if
134  end do
135 
136  npc=0
137  npcm3=0
138 
139 !
140 !.....SEARCH BETWEEN EACH PAIR OF ATOMS
141 !
142 
143  if (stwo) then
144    do jii=1,nb
145      do ii=1,3
146        v1(ii)=xatm(ii,inxat)
147        v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii))
148        vt(ii)=(v1(ii)+v2(ii))/2._dp
149      end do
150      inter=.true.
151      diff=0._dp
152      pom(:)=vt(:)
153      pom(:)=pom(:)-vi(:)
154      diff=vnorm(pom,0)
155      if (diff > maxcpdst) inter=.false.
156      if (inter) then
157        call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0)
158        if (ires==0) then
159          found=.false.
160          if (npc > 0) then
161            do jj=1,npc
162              pom(:)=vt(:)-pc(:,jj)
163              dist=vnorm(pom,0)
164              if (dist < aim_dtset%dpclim) found=.true.
165            end do
166          end if
167          if (.not.found) then
168            pom(:)=vt(:)
169            call bschg1(pom,-1)
170            pcrb(:,npc+1)=pom(:)
171            pom(:)=pom(:)-vi(:)
172            diff=vnorm(pom,0)
173            if (abs(diff) > maxcpdst) found=.true.
174          end if
175          if (.not.found) then
176            npc=npc+1
177            do jj=1,3
178              pc(jj,npc)=vt(jj)
179              evpc(jj,npc)=ev(jj)
180              do kk=1,3
181                zpc(kk,jj,npc)=zz(kk,jj)
182              end do
183            end do
184            i1=ev(1)/abs(ev(1))
185            i2=ev(2)/abs(ev(2))
186            i3=ev(3)/abs(ev(3))
187            icpc(npc)=i1+i2+i3
188            if (icpc(npc)==-3) then
189              npcm3=npcm3+1
190            end if
191            write(std_out,*) 'New critical point found'
192            write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3)
193            write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3)
194            write(std_out,'("AUTOVAL: ",3F16.8)') ev
195            write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') &
196 &           ((zpc(ii,jj,npc),ii=1,3),jj=1,3)
197            call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0)
198            write(22,'(":PC2",3F10.6,3E12.4,I4,2E12.4)') &
199 &           (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
200            write(std_out,'(":PC2",3F10.6,3E12.4,I4,2E12.4)')  &
201 &           (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
202            pom(:)=vt(:)-v1(:)
203            dif1=vnorm(pom,0)
204            pom(:)=vt(:)-v2(:)
205            dif2=vnorm(pom,0)
206            write(std_out,'(":DISPC ",2F12.8)') dif1,dif2
207          end if
208        end if
209      end if
210    end do
211  end if
212 !
213 !.....SEARCH BETWEEN EACH THREE ATOMS
214 !
215  if(sthree) then
216    do jii=1,nb
217      do kjj=jii+1,nb
218        do ii=1,3
219          v1(ii)=xatm(ii,inxat)
220          v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii))
221          v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj))
222          vt(ii)=(v1(ii)+v2(ii)+v3(ii))/3._dp
223        end do
224        inter=.true.
225        pom(:)=vt(:)
226        pom(:)=pom(:)-vi(:)
227        dist=vnorm(pom,0)
228        if (abs(diff)>maxcpdst) then
229          inter=.false.
230          exit
231        end if
232        if (inter) then
233          do jj=1,npc
234            pom(:)=pc(:,jj)-vt(:)
235            diff=vnorm(pom,0)
236            if (diff<aim_dpc0) then
237              inter=.false.
238              exit
239            end if
240          end do
241        end if
242        if (inter) then
243          call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0)
244          if (ires==0) then
245            found=.false.
246            if (npc>0) then
247              do jj=1,npc
248                pom(:)=vt(:)-pc(:,jj)
249                dist=vnorm(pom,0)
250                if (dist<aim_dtset%dpclim) then
251                  found=.true.
252                  exit
253                end if
254              end do
255            end if
256            if (.not.found) then
257              pom(:)=vt(:)
258              call bschg1(pom,-1)
259              pcrb(:,npc+1)=pom(:)
260              pom(:)=pom(:)-vi(:)
261              diff=vnorm(pom,0)
262              if (abs(diff)>maxcpdst) found=.true.
263            end if
264            if (.not.found) then
265              npc=npc+1
266              do jj=1,3
267                pc(jj,npc)=vt(jj)
268                evpc(jj,npc)=ev(jj)
269                do kk=1,3
270                  zpc(kk,jj,npc)=zz(kk,jj)
271                end do
272              end do
273              i1=ev(1)/abs(ev(1))
274              i2=ev(2)/abs(ev(2))
275              i3=ev(3)/abs(ev(3))
276              icpc(npc)=i1+i2+i3
277              if (icpc(npc)==-3) then
278                npcm3=npcm3+1
279              end if
280              write(std_out,*) 'New critical point found'
281              write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3)
282              write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3)
283              write(std_out,'("AUTOVAL: ",3F16.8)') ev
284              write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') &
285 &             ((zpc(ii,jj,npc),ii=1,3),jj=1,3)
286              call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0)
287              write(22,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') &
288 &             (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
289              write(std_out,'(":PC3",3F10.6,3E12.4,I4,2E12.4)') &
290 &             (pc(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
291            end if
292          end if
293        end if
294      end do
295    end do
296  end if
297 
298 !
299 !.....SEARCH BETWEEN EACH FOUR ATOMS
300 !
301  if (sfour) then
302    do jii=1,nb
303      do kjj=jii+1,nb
304        do ll=jii+1,nb
305          do ii=1,3
306            v1(ii)=xatm(ii,inxat)
307            v2(ii)=xatm(ii,ibat(jii))+atp(ii,ipibat(jii))
308            v3(ii)=xatm(ii,ibat(kjj))+atp(ii,ipibat(kjj))
309            v4(ii)=xatm(ii,ibat(ll))+atp(ii,ipibat(ll))
310            vt(ii)=(v1(ii)+v2(ii)+v3(ii)+v4(ii))/4._dp
311          end do
312          inter=.true.
313          pom(:)=vt(:)
314          pom(:)=pom(:)-vi(:)
315          diff=vnorm(pom,0)
316          if (abs(diff)>maxcpdst) then
317            inter=.false.
318            exit
319          end if
320          if (inter) then
321            do jj=1,npc
322              pom(:)=pc(:,jj)-vt(:)
323              diff=vnorm(pom,0)
324              if (diff < aim_dpc0) then
325                inter=.false.
326                exit
327              end if
328            end do
329          end if
330          if (inter) then
331            call critic(aim_dtset,vt,ev,zz,aim_dmaxcs,ires,0)
332            if (ires==0) then
333              found=.false.
334              if (npc>0) then
335                do jj=1,npc
336                  pom(:)=vt(:)-pc(:,jj)
337                  dist=vnorm(pom,0)
338                  if (dist < aim_dtset%dpclim) found=.true.
339                end do
340              end if
341              if (.not.found) then
342                pom(:)=vt(:)
343                pcrb(:,npc+1)=pom(:)
344                pom(:)=pom(:)-vi(:)
345                diff=vnorm(pom,0)
346                if (abs(diff)>maxcpdst) found=.true.
347              end if
348              if (.not.found) then
349                npc=npc+1
350                do jj=1,3
351                  pc(jj,npc)=vt(jj)
352                  evpc(jj,npc)=ev(jj)
353                  do kk=1,3
354                    zpc(kk,jj,npc)=zz(kk,jj)
355                  end do
356                end do
357                i1=ev(1)/abs(ev(1))
358                i2=ev(2)/abs(ev(2))
359                i3=ev(3)/abs(ev(3))
360                icpc(npc)=i1+i2+i3
361                if (icpc(npc)==-3) then
362                  npcm3=npcm3+1
363                end if
364                write(std_out,*) 'New critical point found'
365                write(std_out,'("POS: ",3F16.8)') (pc(ii,npc),ii=1,3)
366                write(std_out,'("POS in base: ",3F16.8)') (pcrb(ii,npc),ii=1,3)
367                write(std_out,'("AUTOVAL: ",3F16.8)') ev
368                write(std_out,'("AUTOVEC: ",/,3F16.8,/,3F16.8,/,3F16.8)') &
369 &               ((zpc(ii,jj,npc),ii=1,3),jj=1,3)
370                call vgh_rho(vt,chg,grho,hrho,rr,iat,ipos,0)
371                write(22,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') &
372 &               (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
373                write(std_out,'(":PC4",3F10.6,3E12.4,I4,2E12.4)') &
374 &               (pc(jj,npc),jj=1,3), (ev(jj),jj=1,3),icpc(npc),ev(1)+ev(2)+ev(3),chg
375              end if
376            end if
377          end if
378        end do
379      end do
380    end do
381  end if
382 
383  write(std_out,*) npc
384 end subroutine critics