hdiff output

r32054/fedor.f 2017-03-08 18:30:08.202158947 +0000 r32053/fedor.f 2017-03-08 18:30:08.462162385 +0000
 11: C   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: C   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: C   GNU General Public License for more details. 13: C   GNU General Public License for more details.
 14: C 14: C
 15: C   You should have received a copy of the GNU General Public License 15: C   You should have received a copy of the GNU General Public License
 16: C   along with this program; if not, write to the Free Software 16: C   along with this program; if not, write to the Free Software
 17: C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: C 18: C
 19:  19: 
 20:            SUBROUTINE ASAR1(R,Ex,Ex1) 20:            SUBROUTINE ASAR1(R,Ex,Ex1)
 21: c          REAL*8  FUNCTION UTN(R,IP) 21: c          REAL*8  FUNCTION UTN(R,I,IP)
 22: C  TERM X0G+ AR2 (AZIZ - JCP 99 (1993) 4518) 22: C  TERM X0G+ AR2 (AZIZ - JCP 99 (1993) 4518)
 23:       implicit DOUBLE PRECISION (a-h,o-z) 23:       implicit DOUBLE PRECISION (a-h,o-z)
 24:       IMPLICIT INTEGER(I-N) 24:       IMPLICIT INTEGER(I-N)
 25:       save 25:       save
 26:       DATA E/143.235d0/,Rm/3.757d0/, A/9.03228328d0/,B/-2.37132823d0/, 26:       DATA E/143.235d0/,Rm/3.757d0/, A/9.03228328d0/,B/-2.37132823d0/,
 27:      * C6,C8,C10/1.09309955d0,.51568309d0,.32521242d0/,AA/8.73933927D4/ 27:      * C6,C8,C10/1.09309955d0,.51568309d0,.32521242d0/,AA/8.73933927D4/
 28:      *,C12,C14/.27818156d0,.31111959d0/,ro/1.107d0/ 28:      *,C12,C14/.27818156d0,.31111959d0/,ro/1.107d0/
 29:      *,NCALL/1/,Ekau/3.1668d-6/,Rau/0.52918d0/ 29:      *,NCALL/1/,Ekau/3.1668d-6/,Rau/0.52918d0/
 30:       IF(NCALL.EQ.1) THEN 30:       IF(NCALL.EQ.1) THEN
 31:       Rm=Rm/Rau 31:       Rm=Rm/Rau
237:         END237:         END
238: 238: 
239: 239: 
240:       subroutine rgnii(n,x,grad,ereal,gradt) !  ,h0,h1,ee,ev,w,natoms))240:       subroutine rgnii(n,x,grad,ereal,gradt) !  ,h0,h1,ee,ev,w,natoms))
241:       implicit DOUBLE PRECISION (a-h,o-z)241:       implicit DOUBLE PRECISION (a-h,o-z)
242:       IMPLICIT INTEGER(I-N)242:       IMPLICIT INTEGER(I-N)
243:       save243:       save
244:         DOUBLE PRECISION x(3,n),h0(9*n*(n-1)/2,9*n*(n-1)/2)244:         DOUBLE PRECISION x(3,n),h0(9*n*(n-1)/2,9*n*(n-1)/2)
245:      &        ,h1(9*n*(n-1)/2,9*n*(n-1)/2)245:      &        ,h1(9*n*(n-1)/2,9*n*(n-1)/2)
246:      &        ,grad(n*3),ee(9*n*(n-1)/2),ev(9*n*(n-1)/2,9*n*(n-1)/2)246:      &        ,grad(n*3),ee(9*n*(n-1)/2),ev(9*n*(n-1)/2,9*n*(n-1)/2)
247: C     &        ,w(9*n*(n-1)/2,1)247:      &        ,w(9*n*(n-1)/2,1)
248:       logical gradt248:       logical gradt
249:         INTEGER INFO, NCALL249:         INTEGER INFO, NCALL
250:         DOUBLE PRECISION TEMPA(9*N)250:         DOUBLE PRECISION TEMPA(9*N)
251:       data iev/1/, ncall/1/, shift/10.0d0/251:       data iev/1/, ncall/1/, shift/10.0d0/
252:       n3=9*n*(n-1)/2252:       n3=9*n*(n-1)/2
253: 253: 
254:         id=1254:         id=1
255:         it=0255:         it=0
256:       call hmatd(x,n,h0,h1,0)256:       call hmatd(x,n,n3,h0,h1,h2,0)
257:       if(ncall.eq.1.or.it.eq.0) then257:       if(ncall.eq.1.or.it.eq.0) then
258:          if(id.eq.0) then258:          if(id.eq.0) then
259: C            call jacob2(h0,ev,w(1,1),ee,w(1,2),n3,iev)259: C            call jacob2(h0,ev,w(1,1),ee,w(1,2),n3,iev)
260:                      else260:                      else
261: C       call tred2(h0,n3,n3,ee,w(1,1),iev)261: C       call tred2(h0,n3,n3,ee,w(1,1),iev)
262: C       call tqli(ee,w(1,1),n3,n3,h0,iev)262: C       call tqli(ee,w(1,1),n3,n3,h0,iev)
263: 263: 
264:         CALL DSYEV('V','U',N3,H0,N3,EE,TEMPA,9*N,INFO)264:         CALL DSYEV('V','U',N3,H0,N3,EE,TEMPA,9*N,INFO)
265:         IF (INFO.NE.0) THEN265:         IF (INFO.NE.0) THEN
266:            PRINT*,'WARNING - INFO=',INFO,' in DSYEV'266:            PRINT*,'WARNING - INFO=',INFO,' in DSYEV'
279:       do i=1,n3279:       do i=1,n3
280:       h0(i,i)=h0(i,i)-shift280:       h0(i,i)=h0(i,i)-shift
281:       enddo281:       enddo
282: C       call iteig(h0,ev(1,n3),ereal,w(1,1),w(1,2),n3)282: C       call iteig(h0,ev(1,n3),ereal,w(1,1),w(1,2),n3)
283:       ereal=ereal+shift283:       ereal=ereal+shift
284:                                 end if284:                                 end if
285: 285: 
286:       if(gradt) then286:       if(gradt) then
287: 287: 
288:       do i=1,n*3288:       do i=1,n*3
289:       call hmatd(x,n,h0,h1,i)289:       call hmatd(x,n,n3,h0,h1,h2,i)
290:       vhvi=0.0d0290:       vhvi=0.0d0
291:       do j=1,n3291:       do j=1,n3
292:       vh=0.0d0292:       vh=0.0d0
293:       do k=1,n3293:       do k=1,n3
294:       vh=vh+ev(k,n3)*h1(k,j)294:       vh=vh+ev(k,n3)*h1(k,j)
295:       enddo295:       enddo
296:       vhvi=vhvi+vh*ev(j,n3)296:       vhvi=vhvi+vh*ev(j,n3)
297:       enddo297:       enddo
298:       grad(i)=vhvi298:       grad(i)=vhvi
299:       enddo299:       enddo
301:               end if301:               end if
302:       return302:       return
303:       end303:       end
304: c ------------------------------------------------------------304: c ------------------------------------------------------------
305: c      real*8 function energy(r0,i,ip)305: c      real*8 function energy(r0,i,ip)
306: c      real*8 r0306: c      real*8 r0
307: c      r=r0307: c      r=r0
308: c      energy=ut(r,i,ip)308: c      energy=ut(r,i,ip)
309: c      return309: c      return
310: c      end310: c      end
311:       subroutine hmatd(x,n,h0,h1,igrad)311:       subroutine hmatd(x,n,n3,h0,h1,h2,igrad)
312:       USE commons312:       USE commons
313:       implicit DOUBLE PRECISION (a-h,o-z)313:       implicit DOUBLE PRECISION (a-h,o-z)
314:       IMPLICIT INTEGER(I-N)314:       IMPLICIT INTEGER(I-N)
315:       INTEGER N315:       INTEGER N, N3
316:       save316:       save
317: C        common/test/itest317: C        common/test/itest
318:       DOUBLE PRECISION x(3,n),h0(9*n*(n-1)/2,9*n*(n-1)/2)318:       DOUBLE PRECISION x(3,n),h0(9*n*(n-1)/2,9*n*(n-1)/2)
319:      1               ,h1(9*n*(n-1)/2,9*n*(n-1)/2)319:      1               ,h1(9*n*(n-1)/2,9*n*(n-1)/2)
320:       INTEGER J0, K, KI, KI0, IJKD, NN9, KJ0, IR, IJK0, IGRAD, I, IJ, IJ0, IJK, J320:       INTEGER J0, K, KI, KI0, IJKD, NN9, KJ0, IR, IJK0, IGRAD, I, IJ, IJ0, IJK, J
321: 321: 
322:       data RauA/0.52918d0/ ,Runit/1.d0/322:       data RauA/0.52918d0/ ,Runit/1.d0/
323: c           isu/1/,isg/2/,ipu/3/,ipg/4/,ix/5/323: c           isu/1/,isg/2/,ipu/3/,ipg/4/,ix/5/
324: 324: 
325: 325: 
338:       enddo338:       enddo
339:       enddo339:       enddo
340: 340: 
341:       ExS=0.0d0341:       ExS=0.0d0
342:       do k=1,n-1342:       do k=1,n-1
343:       do l=k+1,n343:       do l=k+1,n
344:       rkl=344:       rkl=
345:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)345:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
346: cc      Ex=energy(rkl,ix,0)346: cc      Ex=energy(rkl,ix,0)
347:         IF (NEON) THEN347:         IF (NEON) THEN
348:            Ex=utn(rkl/RauA,0)348:            Ex=utn(rkl/RauA,5,0)
349:         ELSE349:         ELSE
350:            CALL ASAR1(rkl,Ex,Ex1)350:            CALL ASAR1(rkl,Ex,Ex1)
351:         ENDIF351:         ENDIF
352:       ExS=ExS+Ex352:       ExS=ExS+Ex
353:       enddo353:       enddo
354:       enddo354:       enddo
355: 355: 
356:       do i=1,n-1356:       do i=1,n-1
357:       do j=i+1,n357:       do j=i+1,n
358: 358: 
359:       ij = (2*n-i)*(i-1)/2 +j-i359:       ij = (2*n-i)*(i-1)/2 +j-i
360:       rij=360:       rij=
361:      1      sqrt((x(1,i)-x(1,j))**2+(x(2,i)-x(2,j))**2+(x(3,i)-x(3,j))**2)361:      1      sqrt((x(1,i)-x(1,j))**2+(x(2,i)-x(2,j))**2+(x(3,i)-x(3,j))**2)
362: cc      Exij=energy(rij,ix,0)362: cc      Exij=energy(rij,ix,0)
363:         IF (NEON) THEN363:         IF (NEON) THEN
364:            Exij=utn(rij/RauA,0)364:            Exij=utn(rij/RauA,5,0)
365:         ELSE365:         ELSE
366:            CALL ASAR1(rij,Exij,Exij1)366:            CALL ASAR1(rij,Exij,Exij1)
367:         ENDIF367:         ENDIF
368: 368: 
369:       do ki=1,3369:       do ki=1,3
370:       do kj=1,3370:       do kj=1,3
371: 371: 
372:       ijk = 9*(ij-1) +3*(ki-1) +kj372:       ijk = 9*(ij-1) +3*(ki-1) +kj
373:       h0(ijk,ijk) = h0(ijk,ijk) +ExS -Exij +Runit/rij373:       h0(ijk,ijk) = h0(ijk,ijk) +ExS -Exij +Runit/rij
374: 374: 
644:                   end if644:                   end if
645: 645: 
646:       ExS1=0.0d0646:       ExS1=0.0d0
647:       do k=1,n-1647:       do k=1,n-1
648:       do l=k+1,n648:       do l=k+1,n
649:       if(l.eq.m.or.k.eq.m) then649:       if(l.eq.m.or.k.eq.m) then
650:       rkl=650:       rkl=
651:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)651:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
652: cc      Ex1= energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl652: cc      Ex1= energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl
653:         IF (NEON) THEN653:         IF (NEON) THEN
654:         Ex1= utn(rkl/RauA,1)/RauA*(x(ir,l)-x(ir,k))/rkl654:         Ex1= utn(rkl/RauA,5,1)/RauA*(x(ir,l)-x(ir,k))/rkl
655:         ELSE655:         ELSE
656:            CALL ASAR1(rkl,Ex,Ex1)656:            CALL ASAR1(rkl,Ex,Ex1)
657:         Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl657:         Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl
658:         ENDIF658:         ENDIF
659:       if(k.eq.m) Ex1=-Ex1659:       if(k.eq.m) Ex1=-Ex1
660:       ExS1=ExS1+Ex1660:       ExS1=ExS1+Ex1
661:                        end if661:                        end if
662:       enddo662:       enddo
663:       enddo663:       enddo
664: 664: 
668:       ij = (2*n-i)*(i-1)/2 +j-i668:       ij = (2*n-i)*(i-1)/2 +j-i
669: 669: 
670:       if(i.eq.m.or.j.eq.m) then670:       if(i.eq.m.or.j.eq.m) then
671: 671: 
672:       rij=672:       rij=
673:      1      sqrt((x(1,i)-x(1,j))**2+(x(2,i)-x(2,j))**2+(x(3,i)-x(3,j))**2)673:      1      sqrt((x(1,i)-x(1,j))**2+(x(2,i)-x(2,j))**2+(x(3,i)-x(3,j))**2)
674:       drijdx=(x(ir,i)-x(ir,j))/rij674:       drijdx=(x(ir,i)-x(ir,j))/rij
675:       if(j.eq.m) drijdx=-drijdx675:       if(j.eq.m) drijdx=-drijdx
676: cc      Exij1= energy(rij,ix,1)*drijdx676: cc      Exij1= energy(rij,ix,1)*drijdx
677:         IF (NEON) THEN677:         IF (NEON) THEN
678:         Exij1= utn(rij/RauA,1)/RauA*drijdx678:         Exij1= utn(rij/RauA,5,1)/RauA*drijdx
679:         ELSE679:         ELSE
680:            CALL ASAR1(rij,Exij,Exij1)680:            CALL ASAR1(rij,Exij,Exij1)
681:         Exij1= Exij1*drijdx681:         Exij1= Exij1*drijdx
682:         ENDIF682:         ENDIF
683: 683: 
684:                            end if684:                            end if
685: 685: 
686:       do ki=1,3686:       do ki=1,3
687:       do kj=1,3687:       do kj=1,3
688: 688: 
1202: 1202: 
1203:       save1203:       save
1204:       data RauA/0.52918d0/1204:       data RauA/0.52918d0/
1205: 1205: 
1206:       ExS=0.0d01206:       ExS=0.0d0
1207:       do k=1,n-11207:       do k=1,n-1
1208:       do l=k+1,n1208:       do l=k+1,n
1209:       rkl=1209:       rkl=
1210:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)1210:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
1211:         IF (NEON) THEN1211:         IF (NEON) THEN
1212:            Ex=utn(rkl/RauA,0)1212:            Ex=utn(rkl/RauA,5,0)
1213:         ELSE1213:         ELSE
1214:            CALL ASAR1(rkl,Ex,Ex1)1214:            CALL ASAR1(rkl,Ex,Ex1)
1215:         ENDIF1215:         ENDIF
1216:       ExS=ExS+Ex1216:       ExS=ExS+Ex
1217:       enddo1217:       enddo
1218:       enddo1218:       enddo
1219:       ereal=ExS1219:       ereal=ExS
1220: 1220: 
1221: 1221: 
1222:       if(gradt) then1222:       if(gradt) then
1225:       do ir=1,31225:       do ir=1,3
1226:       ic=(m-1)*3+ir1226:       ic=(m-1)*3+ir
1227: 1227: 
1228:       ExS1=0.0d01228:       ExS1=0.0d0
1229:       do k=1,n-11229:       do k=1,n-1
1230:       do l=k+1,n1230:       do l=k+1,n
1231:       if(l.eq.m.or.k.eq.m) then1231:       if(l.eq.m.or.k.eq.m) then
1232:       rkl=1232:       rkl=
1233:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)1233:      1      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
1234:         IF (NEON) THEN1234:         IF (NEON) THEN
1235:       Ex1= utn(rkl/RauA,1)/RauA*(x(ir,l)-x(ir,k))/rkl1235:       Ex1= utn(rkl/RauA,5,1)/RauA*(x(ir,l)-x(ir,k))/rkl
1236:       if(k.eq.m) Ex1=-Ex11236:       if(k.eq.m) Ex1=-Ex1
1237:         ELSE1237:         ELSE
1238:            CALL ASAR1(rkl,Ex,Ex1)1238:            CALL ASAR1(rkl,Ex,Ex1)
1239:         Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl1239:         Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl
1240:         if(k.eq.m) Ex1=-Ex11240:         if(k.eq.m) Ex1=-Ex1
1241:         ENDIF1241:         ENDIF
1242:       ExS1=ExS1+Ex11242:       ExS1=ExS1+Ex1
1243:                        end if1243:                        end if
1244:       enddo1244:       enddo
1245:       enddo1245:       enddo
1715: 1715: 
1716:         r0=r1716:         r0=r
1717:         Esu=ut(r0,1,ut1)1717:         Esu=ut(r0,1,ut1)
1718:         Esu1=ut11718:         Esu1=ut1
1719:         Esg=ut(r0,2,ut1)1719:         Esg=ut(r0,2,ut1)
1720:         Esg1=ut11720:         Esg1=ut1
1721:         Epu=ut(r0,3,ut1)1721:         Epu=ut(r0,3,ut1)
1722:         Epu1=ut11722:         Epu1=ut1
1723:         Epg=ut(r0,4,ut1)1723:         Epg=ut(r0,4,ut1)
1724:         Epg1=ut11724:         Epg1=ut1
1725:       Ex =utn(r/RauA,0)1725:       Ex =utn(r/RauA,5,0)
1726:       Ex1=utn(r/RauA,1)/RauA1726:       Ex1=utn(r/RauA,5,1)/RauA
1727:         return1727:         return
1728:         end1728:         end
1729: 1729: 
1730:       DOUBLE PRECISION function ut(rr,i0,ut1)1730:       DOUBLE PRECISION function ut(rr,i0,ut1)
1731:       implicit DOUBLE PRECISION (a-h,o-z)1731:       implicit DOUBLE PRECISION (a-h,o-z)
1732:       IMPLICIT INTEGER(I-N)1732:       IMPLICIT INTEGER(I-N)
1733:       save1733:       save
1734:       common/points/ru,Eu, rg,Eg, rpu,Epu, rpg,Epg, rx,Ex,1734:       common/points/ru,Eu, rg,Eg, rpu,Epu, rpg,Epg, rx,Ex,
1735:      1              coefu,coefg,coefpu,coefpg,coefx, ir1,ir2,ir3,ir4,ir51735:      1              coefu,coefg,coefpu,coefpg,coefx, ir1,ir2,ir3,ir4,ir5
1736:       character(LEN=12) name1,name2,name3,name4,name51736:       character(LEN=12) name1,name2,name3,name4,name5
1796:                     end if1796:                     end if
1797:       if(ind1.eq.2) then1797:       if(ind1.eq.2) then
1798:       ut=ppvalu(ru,coefu,ir1,4,rr,0)1798:       ut=ppvalu(ru,coefu,ir1,4,rr,0)
1799:         ut1=ppvalu(ru,coefu,ir1,4,rr,1) 1799:         ut1=ppvalu(ru,coefu,ir1,4,rr,1) 
1800:                 end if1800:                 end if
1801:                         else1801:                         else
1802:       ut=Eu(ir1) -Easu1802:       ut=Eu(ir1) -Easu
1803:       ut1=0.01803:       ut1=0.0
1804:                         end if1804:                         end if
1805:                     end if1805:                     end if
1806:       if(ind1.eq.1) ut=utn(rr/rauA,0)1806:       if(ind1.eq.1) ut=utn(rr/rauA,i0,0)
1807:       return1807:       return
1808:                  end if1808:                  end if
1809: c1809: c
1810:       if(i0.eq.2) then1810:       if(i0.eq.2) then
1811:       if(ind2.eq.0.or.ind2.eq.2) then1811:       if(ind2.eq.0.or.ind2.eq.2) then
1812:       if(nc2.eq.1) then1812:       if(nc2.eq.1) then
1813:       nc2=21813:       nc2=2
1814:       open(22,file=name2,STATUS='OLD')1814:       open(22,file=name2,STATUS='OLD')
1815:       ir2=01815:       ir2=0
1816:    20 continue1816:    20 continue
1838:                     end if1838:                     end if
1839:       if(ind2.eq.2) then1839:       if(ind2.eq.2) then
1840:       ut=ppvalu(rg,coefg,ir2,4,rr,0)1840:       ut=ppvalu(rg,coefg,ir2,4,rr,0)
1841:         ut1=ppvalu(rg,coefg,ir2,4,rr,1) 1841:         ut1=ppvalu(rg,coefg,ir2,4,rr,1) 
1842:                 end if1842:                 end if
1843:                         else1843:                         else
1844:       ut=Eg(ir2) -Easg1844:       ut=Eg(ir2) -Easg
1845:       ut1=0.01845:       ut1=0.0
1846:                         end if1846:                         end if
1847:                     end if 1847:                     end if 
1848:       if(ind2.eq.1) ut=utn(rr/rauA,0)1848:       if(ind2.eq.1) ut=utn(rr/rauA,i0,0)
1849:       return1849:       return
1850:                  end if1850:                  end if
1851: c1851: c
1852:       if(i0.eq.3) then1852:       if(i0.eq.3) then
1853:       if(ind3.eq.0.or.ind3.eq.2) then1853:       if(ind3.eq.0.or.ind3.eq.2) then
1854:       if(nc3.eq.1) then1854:       if(nc3.eq.1) then
1855:       nc3=21855:       nc3=2
1856:       open(33,file=name3,STATUS='OLD')1856:       open(33,file=name3,STATUS='OLD')
1857:       ir3=01857:       ir3=0
1858:    30 continue1858:    30 continue
1881:                     end if1881:                     end if
1882:       if(ind3.eq.2) then1882:       if(ind3.eq.2) then
1883:       ut=ppvalu(rpu,coefpu,ir3,4,rr,0) 1883:       ut=ppvalu(rpu,coefpu,ir3,4,rr,0) 
1884:         ut1=ppvalu(rpu,coefpu,ir3,4,rr,1) 1884:         ut1=ppvalu(rpu,coefpu,ir3,4,rr,1) 
1885:                 end if1885:                 end if
1886:                         else1886:                         else
1887:       ut=Epu(ir3) -Easpu1887:       ut=Epu(ir3) -Easpu
1888:       ut1=0.01888:       ut1=0.0
1889:                         end if1889:                         end if
1890:                     end if1890:                     end if
1891:       if(ind3.eq.1) ut=utn(rr/rauA,0)1891:       if(ind3.eq.1) ut=utn(rr/rauA,i0,0)
1892:       return1892:       return
1893:                end if1893:                end if
1894: c1894: c
1895:       if(i0.eq.4) then1895:       if(i0.eq.4) then
1896:       if(ind4.eq.0.or.ind4.eq.2) then1896:       if(ind4.eq.0.or.ind4.eq.2) then
1897:       if(nc4.eq.1) then1897:       if(nc4.eq.1) then
1898:       nc4=21898:       nc4=2
1899:       open(44,file=name4,STATUS='OLD')1899:       open(44,file=name4,STATUS='OLD')
1900:       ir4=01900:       ir4=0
1901:    40 continue1901:    40 continue
1924:                     end if1924:                     end if
1925:       if(ind4.eq.2) then1925:       if(ind4.eq.2) then
1926:       ut=ppvalu(rpg,coefpg,ir4,4,rr,0)1926:       ut=ppvalu(rpg,coefpg,ir4,4,rr,0)
1927:         ut1=ppvalu(rpg,coefpg,ir4,4,rr,1)1927:         ut1=ppvalu(rpg,coefpg,ir4,4,rr,1)
1928:                 end if1928:                 end if
1929:                         else1929:                         else
1930:       ut=Epg(ir4) -Easpg1930:       ut=Epg(ir4) -Easpg
1931:       ut1=0.01931:       ut1=0.0
1932:                         end if1932:                         end if
1933:                     end if 1933:                     end if 
1934:       if(ind4.eq.1) ut=utn(rr/rauA,0)1934:       if(ind4.eq.1) ut=utn(rr/rauA,i0,0)
1935:       return1935:       return
1936:                  end if1936:                  end if
1937: c1937: c
1938:       if(i0.eq.5) then1938:       if(i0.eq.5) then
1939:       if(ind5.eq.0.or.ind5.eq.2) then1939:       if(ind5.eq.0.or.ind5.eq.2) then
1940:       if(nc5.eq.1) then1940:       if(nc5.eq.1) then
1941:       nc5=21941:       nc5=2
1942:       open(55,file=name5,STATUS='OLD')1942:       open(55,file=name5,STATUS='OLD')
1943:       ir5=01943:       ir5=0
1944:    50 continue1944:    50 continue
1966:                     end if1966:                     end if
1967:       if(ind5.eq.2) then1967:       if(ind5.eq.2) then
1968:       ut=ppvalu(rx,coefx,ir5,4,rr,0)1968:       ut=ppvalu(rx,coefx,ir5,4,rr,0)
1969:       ut1=ppvalu(rx,coefx,ir5,4,rr,1)1969:       ut1=ppvalu(rx,coefx,ir5,4,rr,1)
1970:                 end if1970:                 end if
1971:                         else 1971:                         else 
1972:       ut=Ex(ir5) -Easx1972:       ut=Ex(ir5) -Easx
1973:       ut1=0.01973:       ut1=0.0
1974:                         end if1974:                         end if
1975:                     end if1975:                     end if
1976:       if(ind5.eq.1) ut=utn(rr/rauA,0)1976:       if(ind5.eq.1) ut=utn(rr/rauA,i0,0)
1977:       return1977:       return
1978:                  end if1978:                  end if
1979:       return 1979:       return 
1980:       end1980:       end
1981: 1981: 
1982: 1982: 
1983:       subroutine rgni(n,x,grad,ereal,gradt) !  ,h0,h1,ee,ev,w,natoms)1983:       subroutine rgni(n,x,grad,ereal,gradt) !  ,h0,h1,ee,ev,w,natoms)
1984:         implicit DOUBLE PRECISION (a-h,o-z)1984:         implicit DOUBLE PRECISION (a-h,o-z)
1985:       IMPLICIT INTEGER(I-N)1985:       IMPLICIT INTEGER(I-N)
1986:         INTEGER INFO, NEMAX1986:         INTEGER INFO, NEMAX
1987:         DOUBLE PRECISION TEMPA(9*N)1987:         DOUBLE PRECISION TEMPA(9*N)
1988: 1988: 
1989:       parameter (nemax=100)1989:       parameter (nemax=100)
1990:         integer iw(nemax,3), NCALL, N, I, ID, IEV, IT, J, K, N3, NE1990:         integer iw(nemax,3),  IWORK(3*5*N), IFAIL(3*N), NCALL, N, I, ID, IEV, IT, J, K, N3
1991: C       INTEGER IWORK(3*5*N), IFAIL(3*N) 
1992: 1991: 
1993:         DOUBLE PRECISION x(3,n),h0(n*3,n*3)1992:         DOUBLE PRECISION x(3,n),h0(n*3,n*3)
1994:      &        ,h1(n*3,n*3)1993:      &        ,h1(n*3,n*3)
1995:      &        ,grad(n*3),ee(n*3),ev(n*3,n*3)1994:      &        ,grad(n*3),ee(n*3),ev(n*3,n*3),w(n*3,n*3),WORK(3*8*N)
1996: C     &       ,w(n*3,n*3),WORK(3*8*N) 
1997:       logical gradt1995:       logical gradt
1998:       data iev/1/, ncall/1/, shift/10.0d0/1996:       data iev/1/, ncall/1/, shift/10.0d0/
1999:      &  , nstep/1000/,iflag/-2/1997:      &  , Elow/-1.d3/,eps/1.d-6/,nstep/1000/,message/6/,nwork/777/
2000: C     &  , Elow/-1.d3/,eps/1.d-6/,message/6/,nwork/777/1998:      &  , iflag/-2/ 
2001:       save1999:       save
2002: 2000: 
2003:       ne=0 
2004:       n3=n*32001:       n3=n*3
2005: cc      if(gradt) igrad=12002: cc      if(gradt) igrad=1
2006: 2003: 
2007:         id=12004:         id=1
2008:         it=02005:         it=0
2009: c       it=12006: c       it=1
2010: C       it=22007: C       it=2
2011:       call hmat(x,n,n3,h0,h1,0)2008:       call hmat(x,n,n3,h0,h1,h2,0)
2012:       if(ncall.eq.1.or.it.eq.0) then2009:       if(ncall.eq.1.or.it.eq.0) then
2013:          if(id.eq.0) then2010:          if(id.eq.0) then
2014: C            call jacob2(h0,ev,w(1,1),ee,w(1,2),n3,iev)2011: C            call jacob2(h0,ev,w(1,1),ee,w(1,2),n3,iev)
2015:                      else2012:                      else
2016: C       call tred2(h0,n3,n3,ee,w(1,1),iev)2013: C       call tred2(h0,n3,n3,ee,w(1,1),iev)
2017: C       call tqli(ee,w(1,1),n3,n3,h0,iev)2014: C       call tqli(ee,w(1,1),n3,n3,h0,iev)
2018: 2015: 
2019:         CALL DSYEV('V','U',N3,H0,N3,EE,TEMPA,9*N,INFO)2016:         CALL DSYEV('V','U',N3,H0,N3,EE,TEMPA,9*N,INFO)
2020:         IF (INFO.NE.0) THEN2017:         IF (INFO.NE.0) THEN
2021:            PRINT*,'WARNING - INFO=',INFO,' in DSYEV'2018:            PRINT*,'WARNING - INFO=',INFO,' in DSYEV'
2047:       if(it.eq.1) then2044:       if(it.eq.1) then
2048:       do i=1,n32045:       do i=1,n3
2049:       h0(i,i)=h0(i,i)-shift2046:       h0(i,i)=h0(i,i)-shift
2050:       enddo2047:       enddo
2051: C        call iteig(h0,ev(1,n3),ereal,w(1,1),w(1,2),n3)2048: C        call iteig(h0,ev(1,n3),ereal,w(1,1),w(1,2),n3)
2052:       ereal=ereal+shift2049:       ereal=ereal+shift
2053:                 end if2050:                 end if
2054: 2051: 
2055:       if(it.eq.2) then2052:       if(it.eq.2) then
2056: 2053: 
2057: C        Ehigh=dmax1(0.d0,2.d0*ereal)2054:         Ehigh=dmax1(0.d0,2.d0*ereal)
2058:       do istep=1,nstep2055:       do istep=1,nstep
2059: c      call itlane(n3,Elow,Ehigh,eps,n,n3,nstep,message,nwork2056: c      call itlane(n3,Elow,Ehigh,eps,n,n3,nstep,message,nwork
2060: c     &      ,iflag,ev,ev(1,n3),ee,iw,ne,w,w(1,2),iw(1,3),h1,h1(1,2))2057: c     &      ,iflag,ev,ev(1,n3),ee,iw,ne,w,w(1,2),iw(1,3),h1,h1(1,2))
2061:       if(iflag.eq.0) go to 12058:       if(iflag.eq.0) go to 1
2062:       if(iflag.eq.1) then2059:       if(iflag.eq.1) then
2063:       do i=1,n32060:       do i=1,n3
2064:       hvi=0.d02061:       hvi=0.d0
2065:       do j=1,n32062:       do j=1,n3
2066:       hvi=hvi+h0(i,j)*ev(j,n3)2063:       hvi=hvi+h0(i,j)*ev(j,n3)
2067:       enddo2064:       enddo
2089: C        stop2086: C        stop
2090: C     end if2087: C     end if
2091: 2088: 
2092:                 end if2089:                 end if
2093: 2090: 
2094:                                 end if2091:                                 end if
2095: 2092: 
2096:       if(gradt) then2093:       if(gradt) then
2097: 2094: 
2098:       do i=1,n32095:       do i=1,n3
2099:       call hmat(x,n,n3,h0,h1,i)2096:       call hmat(x,n,n3,h0,h1,h2,i)
2100:       vhvi=0.0d02097:       vhvi=0.0d0
2101:       do j=1,n32098:       do j=1,n3
2102:       vh=0.0d02099:       vh=0.0d0
2103:       do k=1,n32100:       do k=1,n3
2104: cc      vh=vh+ev(k,n3)*h1(k,j,i)2101: cc      vh=vh+ev(k,n3)*h1(k,j,i)
2105:       vh=vh+ev(k,n3)*h1(k,j)2102:       vh=vh+ev(k,n3)*h1(k,j)
2106:       enddo2103:       enddo
2107:       vhvi=vhvi+vh*ev(j,n3)2104:       vhvi=vhvi+vh*ev(j,n3)
2108:       enddo2105:       enddo
2109:       grad(i)=vhvi2106:       grad(i)=vhvi
2139:       do 16 k=1,n2136:       do 16 k=1,n
2140:       evnn=ev(k,nn)2137:       evnn=ev(k,nn)
2141:       ev(k,nn)=ev(k,ind)2138:       ev(k,nn)=ev(k,ind)
2142:       ev(k,ind)=evnn2139:       ev(k,ind)=evnn
2143:    16 continue2140:    16 continue
2144:                    end if2141:                    end if
2145:    17 CONTINUE2142:    17 CONTINUE
2146:       RETURN2143:       RETURN
2147:       END2144:       END
2148: 2145: 
2149:       subroutine hmat(x,n,n3,h0,h1,igrad)2146:       subroutine hmat(x,n,n3,h0,h1,h2,igrad)
2150:       USE commons2147:       USE commons
2151:       implicit DOUBLE PRECISION (a-h,o-z)2148:       implicit DOUBLE PRECISION (a-h,o-z)
2152:       IMPLICIT INTEGER(I-N)2149:       IMPLICIT INTEGER(I-N)
2153: C       common/test/itest2150: C       common/test/itest
2154: 2151: 
2155:       DOUBLE PRECISION x(3,n),h0(n*3,n*3)2152:       DOUBLE PRECISION x(3,n),h0(n*3,n*3)
2156:      &               ,h1(n*3,n*3)2153:      &               ,h1(n*3,n*3)
2157: 2154: 
2158:       save2155:       save
2159:       data RauA/0.52918d0/2156:       data RauA/0.52918d0/
2167:       enddo2164:       enddo
2168:       enddo2165:       enddo
2169: 2166: 
2170:       ExS=0.0d02167:       ExS=0.0d0
2171:       do k=1,n-12168:       do k=1,n-1
2172:       do l=k+1,n2169:       do l=k+1,n
2173:       rkl=2170:       rkl=
2174:      &      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)2171:      &      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
2175: cc      Ex=energy(rkl,ix,0)2172: cc      Ex=energy(rkl,ix,0)
2176:         IF (NEON) THEN2173:         IF (NEON) THEN
2177:            Ex=utn(rkl/RauA,0)2174:            Ex=utn(rkl/RauA,5,0)
2178:         ELSE2175:         ELSE
2179:            CALL ASAR1(rkl,Ex,Ex1)2176:            CALL ASAR1(rkl,Ex,Ex1)
2180:         ENDIF2177:         ENDIF
2181:       ExS=ExS+Ex2178:       ExS=ExS+Ex
2182:       enddo2179:       enddo
2183:       enddo2180:       enddo
2184: 2181: 
2185:       do i=1,n3-12182:       do i=1,n3-1
2186:       h0(i,i)=h0(i,i)+ExS2183:       h0(i,i)=h0(i,i)+ExS
2187:       do j=i+1,n32184:       do j=i+1,n3
2365: 2362: 
2366:       ExS1=0.0d02363:       ExS1=0.0d0
2367:       do k=1,n-12364:       do k=1,n-1
2368:       do l=k+1,n2365:       do l=k+1,n
2369:       if(l.eq.m.or.k.eq.m) then2366:       if(l.eq.m.or.k.eq.m) then
2370:       rkl=2367:       rkl=
2371:      &      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)2368:      &      sqrt((x(1,l)-x(1,k))**2+(x(2,l)-x(2,k))**2+(x(3,l)-x(3,k))**2)
2372: cc      if(l.eq.m) Ex1= energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl2369: cc      if(l.eq.m) Ex1= energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl
2373: cc      if(k.eq.m) Ex1=-energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl2370: cc      if(k.eq.m) Ex1=-energy(rkl,ix,1)*(x(ir,l)-x(ir,k))/rkl
2374:         IF (NEON) THEN2371:         IF (NEON) THEN
2375:       if(l.eq.m) Ex1= utn(rkl/RauA,1)/RauA*(x(ir,l)-x(ir,k))/rkl2372:       if(l.eq.m) Ex1= utn(rkl/RauA,5,1)/RauA*(x(ir,l)-x(ir,k))/rkl
2376:       if(k.eq.m) Ex1=-utn(rkl/RauA,1)/RauA*(x(ir,l)-x(ir,k))/rkl2373:       if(k.eq.m) Ex1=-utn(rkl/RauA,5,1)/RauA*(x(ir,l)-x(ir,k))/rkl
2377:         ELSE2374:         ELSE
2378:            CALL ASAR1(rkl,Ex,Ex1)2375:            CALL ASAR1(rkl,Ex,Ex1)
2379:         if(l.eq.m) Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl2376:         if(l.eq.m) Ex1= Ex1*(x(ir,l)-x(ir,k))/rkl
2380:         if(k.eq.m) Ex1=-Ex1*(x(ir,l)-x(ir,k))/rkl2377:         if(k.eq.m) Ex1=-Ex1*(x(ir,l)-x(ir,k))/rkl
2381:         ENDIF2378:         ENDIF
2382:       ExS1=ExS1+Ex12379:       ExS1=ExS1+Ex1
2383:                        end if2380:                        end if
2384:       enddo2381:       enddo
2385:       enddo2382:       enddo
2386: 2383: 
2648: cc      enddo2645: cc      enddo
2649: cc      enddo2646: cc      enddo
2650: 2647: 
2651:         IF(DIPOLE) CALL dipoles(x,n,h0,h1,igrad)2648:         IF(DIPOLE) CALL dipoles(x,n,h0,h1,igrad)
2652: 2649: 
2653:    2      continue2650:    2      continue
2654: 2651: 
2655:       return2652:       return
2656:       end2653:       end
2657: 2654: 
2658: cc           real*4 function utn(r,ip)2655: cc           real*4 function utn(r,i,ip)
2659: cc           utn=0.0d02656: cc           utn=0.0d0
2660: cc           if(i.ne.5) return2657: cc           if(i.ne.5) return
2661: cc           utn=asne1(r,i,ip)2658: cc           utn=asne1(r,i,ip)
2662: cc           return2659: cc           return
2663: cc           end2660: cc           end
2664: cc                 FUNCTION ASNE1(R,I,IP)2661: cc                 FUNCTION ASNE1(R,I,IP)
2665:                  DOUBLE PRECISION FUNCTION UTN(R,IP)2662:                  DOUBLE PRECISION FUNCTION UTN(R,I,IP)
2666:            implicit DOUBLE PRECISION (a-h,o-z)2663:            implicit DOUBLE PRECISION (a-h,o-z)
2667:       IMPLICIT INTEGER(I-N)2664:       IMPLICIT INTEGER(I-N)
2668:            save2665:            save
2669:              COMMON/DIPOL/ adip,Runit,ReX                                            2666:              COMMON/DIPOL/ adip,Runit,ReX                                            
2670:            DATA E/1.34d-4/,RM/5.841d0/, A/13.86434671d0/,B/-.12993822d0/, 2667:            DATA E/1.34d-4/,RM/5.841d0/, A/13.86434671d0/,B/-.12993822d0/, 
2671:      & D/1.36d0/,2668:      & D/1.36d0/,
2672:      * C6,C8,C10/1.21317545d0,.53222749d0,.24570703d0/, AA/8.9571795d5/2669:      * C6,C8,C10/1.21317545d0,.53222749d0,.24570703d0/, AA/8.9571795d5/
2673:       ReX=RM2670:       ReX=RM
2674:            X=R/RM2671:            X=R/RM
2675:            F=1.d02672:            F=1.d0
2676:            UTN = 0.D0 
2677:            IF(X.LT.D) F=EXP(-(D/X-1.d0)**2)2673:            IF(X.LT.D) F=EXP(-(D/X-1.d0)**2)
2678: 2674: 
2679:            X2=X*X2675:            X2=X*X
2680:            X4=X2*X22676:            X4=X2*X2
2681:            X6=X4*X22677:            X6=X4*X2
2682:            X8=X4*X42678:            X8=X4*X4
2683:            X10=X6*X42679:            X10=X6*X4
2684:            if(ip.eq.0) then2680:            if(ip.eq.0) then
2685:            V=AA*EXP((-A+B*X)*X)-F*(C6/X6+C8/X8+C10/X10)2681:            V=AA*EXP((-A+B*X)*X)-F*(C6/X6+C8/X8+C10/X10)
2686: cc           ASNE1=V*E2682: cc           ASNE1=V*E
2761:       enddo2757:       enddo
2762: 2758: 
2763:       return2759:       return
2764:       end2760:       end
2765: 2761: 
2766: 2762: 
2767:       subroutine ehecl2(r,c,e,e1,e1c)2763:       subroutine ehecl2(r,c,e,e1,e1c)
2768:         implicit DOUBLE PRECISION (a-h,o-z)2764:         implicit DOUBLE PRECISION (a-h,o-z)
2769:       IMPLICIT INTEGER(I-N)2765:       IMPLICIT INTEGER(I-N)
2770:       common/init/Es,Ep,Epb2766:       common/init/Es,Ep,Epb
2771:       data hre/0.995D0/2767:       data hre/0.995D0/, nc/1/
2772: C     &, nc/1/ 
2773:       save2768:       save
2774: 2769: 
2775:       Es = hecl2ut0(r,1,Es1)2770:       Es = hecl2ut0(r,1,Es1)
2776:       Ep = hecl2ut0(r,2,Ep1)2771:       Ep = hecl2ut0(r,2,Ep1)
2777:         Epb= hecl2ut0(r,3,Epb1)2772:         Epb= hecl2ut0(r,3,Epb1)
2778: 2773: 
2779:         ratio=hre/r2774:         ratio=hre/r
2780:         A= 1.d0 + ratio*ratio2775:         A= 1.d0 + ratio*ratio
2781:         B= 2.*ratio2776:         B= 2.*ratio
2782:         A1= -2.d0*(A-1.d0)/r2777:         A1= -2.d0*(A-1.d0)/r


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0