source/unres/src-HCD-5D/dfa.F

   1       subroutine init_dfa_vars
   2
   3       include 'DIMENSIONS'
   4       include 'COMMON.INTERACT'
   5       include 'COMMON.DFA'
   6
   7       integer ii
   8
   9 C     Number of restraints
  10       idisnum = 0
  11       iphinum = 0
  12       ithenum = 0
  13       ineinum = 0
  14
  15       idislis = 0
  16       iphilis = 0
  17       ithelis = 0
  18       ineilis = 0
  19       jneilis = 0
  20       jneinum = 0
  21       kshell  = 0
  22       fnei    = 0
  23 C     For beta
  24       nca     = 0
  25       icaidx  = 0
  26
  27 C     real variables
  28 CC    WEIGHTS for each min
  29       sccdist = 0.0d0
  30       fdist   = 0.0d0
  31       sccphi  = 0.0d0
  32       sccthe  = 0.0d0
  33       sccnei  = 0.0d0
  34       fphi1   = 0.0d0
  35       fphi2   = 0.0d0
  36       fthe1   = 0.0d0
  37       fthe2   = 0.0d0
  38 C     energies
  39       edfatot = 0.0d0
  40       edfadis = 0.0d0
  41       edfaphi = 0.0d0
  42       edfathe = 0.0d0
  43       edfanei = 0.0d0
  44       edfabet = 0.0d0
  45 C     weights for each E term
  46 C     these should be identical with
  47       dis_inc = 0.0d0
  48       phi_inc = 0.0d0
  49       the_inc = 0.0d0
  50       nei_inc = 0.0d0
  51       beta_inc = 0.0d0
  52       wshet   = 0.0d0
  53 C     precalculate exp table!
  54 c      dfaexp  = 0.0d0
  55 c      do ii = 1, 15001
  56 c         dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
  57 c      end do
  58
  59       ishiftca=nnt-1
  60       ilastca=nct
  61
  62       print *,'ishiftca=',ishiftca,'ilastca=',ilastca
  63
  64       return
  65       end
  66
  67
  68       subroutine read_dfa_info
  69 C
  70 C     read fragment informations
  71 C
  72       implicit real*8 (a-h,o-z)
  73       include 'DIMENSIONS'
  74 #ifdef MPI
  75       include "mpif.h"
  76       include 'COMMON.SETUP'
  77       integer ierror
  78 #endif
  79       include 'COMMON.IOUNITS'
  80       include 'COMMON.CHAIN'
  81       include 'COMMON.DFA'
  82       include 'COMMON.FFIELD'
  83       include 'COMMON.CONTROL'
  84
  85
  86 C     NOTE THAT FILENAMES are FIXED, CURRENTLY!!
  87 C     THIS SHOULD BE MODIFIED!!
  88
  89       character*320 buffer
  90       integer iodfa
  91       parameter(iodfa=89)
  92
  93       integer i, j, nval
  94       integer ica1, ica2,ica3,ica4,ica5
  95       integer ishell, inca, itmp,iitmp
  96       double precision wtmp
  97       logical lprn /.false./
  98 C
  99 C     READ DISTANCE
 100 C
 101       open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
 102       goto 34
 103  33   write(iout,'(a)') 'Error opening dist_dfa.dat file'
 104       print *,'Error opening dist_dfa.dat file'
 105       stop
 106  34   continue
 107       write(iout,'(a)') 'dist_dfa.dat is opened!'
 108 C     read title
 109       read(iodfa, '(a)') buffer
 110 C     read number of restraints
 111       read(iodfa, *) IDFADIS
 112       read(iodfa, *) dis_inc
 113       do i=1, idfadis
 114          read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
 115
 116          idisnum(i)=nval
 117          idislis(1,i)=ica1
 118          idislis(2,i)=ica2
 119
 120          do j=1, nval
 121             read(iodfa,*) tmp
 122             fdist(i,j) = tmp
 123          enddo
 124
 125          do j=1, nval
 126             read(iodfa,*) tmp
 127             sccdist(i,j) = tmp
 128          enddo
 129
 130       enddo
 131       close(iodfa)
 132
 133 C     READ ANGLE RESTRAINTS
 134 C     PHI RESTRAINTS
 135       open(iodfa, file='phi_dfa.dat',status='old',err=35)
 136       goto 36
 137  35   write(iout,'(a)') 'Error opening phi_dfa.dat file'
 138       print *, 'Error opening phi_dfa.dat file'
 139       stop
 140
 141  36   continue
 142       write(iout,'(a)') 'phi_dfa.dat is opened!'
 143
 144 C     READ TITLE
 145       read(iodfa, '(a)') buffer
 146 C     READ NUMBER OF RESTRAINTS
 147       READ(iodfa, *) IDFAPHI
 148       read(iodfa,*) phi_inc
 149       do i=1, idfaphi
 150          read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
 151
 152          iphinum(i)=nval
 153
 154          iphilis(1,i)=ica1
 155          iphilis(2,i)=ica2
 156          iphilis(3,i)=ica3
 157          iphilis(4,i)=ica4
 158          iphilis(5,i)=ica5
 159
 160          do j=1, nval
 161             read(iodfa,*) tmp1,tmp2
 162             fphi1(i,j) = tmp1
 163             fphi2(i,j) = tmp2
 164          enddo
 165
 166          do j=1, nval
 167             read(iodfa,*) tmp
 168             sccphi(i,j) = tmp
 169          enddo
 170
 171       enddo
 172       close(iodfa)
 173
 174 C     THETA RESTRAINTS
 175       open(iodfa, file='theta_dfa.dat',status='old',err=41)
 176       goto 42
 177  41   write(iout,'(a)') 'Error opening theta_dfa.dat file'
 178       print *,'Error opening theta_dfa.dat file'
 179       stop
 180  42   continue
 181       write(iout,'(a)') 'theta_dfa.dat is opened!'
 182 C     READ TITLE
 183       read(iodfa, '(a)') buffer
 184 C     READ NUMBER OF RESTRAINTS
 185       READ(iodfa, *) IDFATHE
 186       read(iodfa,*) the_inc
 187
 188       do i=1, idfathe
 189          read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
 190
 191          ithenum(i)=nval
 192
 193          ithelis(1,i)=ica1
 194          ithelis(2,i)=ica2
 195          ithelis(3,i)=ica3
 196          ithelis(4,i)=ica4
 197          ithelis(5,i)=ica5
 198
 199          do j=1, nval
 200             read(iodfa,*) tmp1,tmp2
 201             fthe1(i,j) = tmp1
 202             fthe2(i,j) = tmp2
 203          enddo
 204
 205          do j=1, nval
 206             read(iodfa,*) tmp
 207             sccthe(i,j) = tmp
 208          enddo
 209
 210       enddo
 211       close(iodfa)
 212 C     END of READING ANGLE RESTRAINT!
 213
 214 C     NUMBER OF NEIGHBOR CAs
 215       open(iodfa,file='nei_dfa.dat',status='old',err=37)
 216       goto 38
 217  37   write(iout,'(a)') 'Error opening nei_dfa.dat file'
 218       print *,'Error opening nei_dfa.dat file'
 219       stop
 220  38   continue
 221       write(iout,'(a)') 'nei_dfa.dat is opened!'
 222 C     READ TITLE
 223       read(iodfa, '(a)') buffer
 224 C     READ NUMBER OF RESTRAINTS
 225       READ(iodfa, *) idfanei
 226       read(iodfa,*) nei_inc
 227
 228       do i=1, idfanei
 229          read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
 230
 231          ineilis(i)=ica1
 232          kshell(i)=ishell
 233          ineinum(i)=nval
 234
 235          do j=1, nval
 236             read(iodfa,*) inca
 237             fnei(i,j) = inca
 238 C            write(*,*) 'READ NEI:',i,j,fnei(i,j)
 239          enddo
 240
 241          do j=1, nval
 242             read(iodfa,*) tmp
 243             sccnei(i,j) = tmp
 244          enddo
 245
 246       enddo
 247       close(iodfa)
 248 C     END OF NEIGHBORING CA
 249
 250 C     init parallel
 251 C     BETA is not parallel !
 252 #ifdef MPI
 253       if (wdfa_beta.ne.0.0 .and. nfgtasks.gt.1) then
 254             write (iout,*) "ERRROR dfa_beta works only for FGPROCS=1"
 255             print *,"ERRROR dfa_beta works only for FGPROCS=1"
 256             stop
 257       endif
 258       call int_bounds(idfadis,idfadis_start,idfadis_end)
 259       call int_bounds(idfaphi,idfaphi_start,idfaphi_end)
 260       call int_bounds(idfathe,idfathe_start,idfathe_end)
 261       call int_bounds(idfanei,idfanei_start,idfanei_end)
 262       if (lprn) write (*,*) "Processor",MyRank," DFA MPI ",
 263      & "idfadis ",idfadis,idfadis_start,idfadis_end,
 264      & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end,
 265      & "idfathe ",idfathe,idfathe_start,idfathe_end,
 266      & "idfanei ",idfanei,idfanei_start,idfanei_end
 267       if (lprn)
 268      &  write (iout,*) "DFA MPI ",
 269      & "idfadis ",idfadis,idfadis_start,idfadis_end,
 270      & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end,
 271      & "idfathe ",idfathe,idfathe_start,idfathe_end,
 272      & "idfanei ",idfanei,idfanei_start,idfanei_end
 273       do i=0,max_fg_procs-1
 274         idfadis_start_all(j)=0
 275         idfadis_end_all(j)=0
 276         idfaphi_start_all(j)=0
 277         idfaphi_end_all(j)=0
 278         idfathe_start_all(j)=0
 279         idfathe_end_all(j)=0
 280         idfanei_start_all(j)=0
 281         idfanei_end_all(j)=0
 282       enddo
 283       if (nfgtasks.gt.1) then
 284       call MPI_Allgather(idfadis_start,1,MPI_INTEGER,
 285      &  idfadis_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 286       call MPI_Allgather(idfadis_end,1,MPI_INTEGER,
 287      &  idfadis_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 288       call MPI_Allgather(idfaphi_start,1,MPI_INTEGER,
 289      &  idfaphi_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 290       call MPI_Allgather(idfaphi_end,1,MPI_INTEGER,
 291      &  idfaphi_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 292       call MPI_Allgather(idfathe_start,1,MPI_INTEGER,
 293      &  idfathe_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 294       call MPI_Allgather(idfathe_end,1,MPI_INTEGER,
 295      &  idfathe_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 296       call MPI_Allgather(idfanei_start,1,MPI_INTEGER,
 297      &  idfanei_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 298       call MPI_Allgather(idfanei_end,1,MPI_INTEGER,
 299      &  idfanei_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
 300       if (me.eq.0 .or. out1file) then
 301         write (iout,*) "Partitioning of DFA work"
 302         write (iout,'(5a10)') 'Rank','DFA_dis','DFA_phi','DFA_the',
 303      &   'DFA_nei'
 304         do i=0,nfgtasks-1
 305           write (iout,'(i10,8i5)') i,idfadis_start_all(i),
 306      &     idfadis_end_all(i),idfaphi_start_all(i),
 307      &     idfaphi_end_all(i),idfathe_start_all(i),
 308      &     idfathe_end_all(i),idfanei_start_all(i),
 309      &     idfanei_end_all(i)
 310         enddo
 311       endif
 312       endif
 313 #else
 314       idfadis_start=1
 315       idfadis_end=idfadis
 316       idfaphi_start=1
 317       idfaphi_end=idfaphi
 318       idfathe_start=1
 319       idfathe_end=idfathe
 320       idfanei_start=1
 321       idfanei_end=idfanei
 322 #endif
 323
 324
 325 C     READ BETA RESTRAINT
 326       if (wdfa_beta.eq.0.0) return
 327       open(iodfa, file='beta_dfa.dat',status='old',err=39)
 328       goto 40
 329  39   write(iout,'(a)') 'Error opening beta_dfa.dat file'
 330       print *,'Error opening beta_dfa.dat file'
 331       stop
 332  40   continue
 333       write(iout,'(a)') 'beta_dfa.dat is opened!'
 334
 335       read(iodfa,'(a)') buffer
 336       read(iodfa,*) itmp
 337       read(iodfa,*) beta_inc
 338
 339       do i=1,itmp
 340          read(iodfa,*) ica1, iitmp
 341          do j=1,itmp
 342             read(iodfa,*) wtmp
 343             wshet(i,j) =  wtmp
 344 c            write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
 345          enddo
 346       enddo
 347
 348       close(iodfa)
 349 C     END OF BETA RESTRAINT
 350
 351
 352       return
 353       END
 354
 355       subroutine edfad(edfadis)
 356
 357       implicit real*8 (a-h,o-z)
 358       include 'DIMENSIONS'
 359       include 'COMMON.CHAIN'
 360       include 'COMMON.DERIV'
 361       include 'COMMON.DFA'
 362
 363       double precision edfadis
 364       integer i, iatm1, iatm2,idiff
 365       double precision ckk, sckk,dist,texp
 366       double precision jix,jiy,jiz,ep,fp,scc
 367
 368       edfadis=0
 369       gdfad=0.0d0
 370
 371       do i=idfadis_start,idfadis_end
 372
 373          iatm1=idislis(1,i)+ishiftca
 374          iatm2=idislis(2,i)+ishiftca
 375          idiff = abs(iatm1-iatm2)
 376
 377          JIX=c(1,iatm2)-c(1,iatm1)
 378          JIY=c(2,iatm2)-c(2,iatm1)
 379          JIZ=c(3,iatm2)-c(3,iatm1)
 380          DIST=SQRT(JIX*JIX+JIY*JIY+JIZ*JIZ)
 381
 382          ckk=ck(idiff)
 383          sckk=sck(idiff)
 384
 385          scc = 0.0d0
 386          ep = 0.0d0
 387          fp = 0.0d0
 388
 389          do j=1,idisnum(i)
 390
 391             dd = dist-fdist(i,j)
 392             dtmp = dd*dd/ckk
 393             if (dtmp.ge.15.0d0) then
 394                texp = 0.0d0
 395             else
 396 c               texp = dfaexp( idint(dtmp*1000)+1 )/sckk
 397                 texp = exp(-dtmp)/sckk
 398             endif
 399
 400             ep=ep+sccdist(i,j)*texp
 401             fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
 402             scc=scc+sccdist(i,j)
 403 C            write(*,'(2i8,6f12.5)') i, j, dist,
 404 C     &           fdist(i,j), ep, fp, sccdist(i,j), scc
 405
 406          enddo
 407
 408          ep = -ep/scc
 409          fp = fp/scc
 410
 411
 412 c         IF(ABS(EP).lt.1.0d-20)THEN
 413 c            EP=0.0D0
 414 c         ENDIF
 415 c         IF (ABS(FP).lt.1.0d-20) THEN
 416 c            FP=0.0D0
 417 c         ENDIF
 418
 419          edfadis=edfadis+ep*dis_inc*wwdist
 420
 421          gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
 422          gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
 423          gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
 424
 425          gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
 426          gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
 427          gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
 428
 429       enddo
 430
 431       return
 432       end
 433
 434       subroutine edfat(edfator)
 435 C     DFA torsion angle
 436       implicit real*8 (a-h,o-z)
 437       include 'DIMENSIONS'
 438       include 'COMMON.CHAIN'
 439       include 'COMMON.DERIV'
 440       include 'COMMON.DFA'
 441
 442       integer i,j,ii,iii
 443       integer iatom(5)
 444       double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
 445       double precision cwidth, cwidth2
 446       PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
 447       parameter (TENM20=1.0d-20)
 448
 449       edfator= 0.0d0
 450       enephi = 0.0d0
 451       enethe = 0.0d0
 452       gdfat(:,:) = 0.0d0
 453
 454 C     START OF PHI ANGLE
 455       do i=idfaphi_start,idfaphi_end
 456
 457          aphi = 0.0d0
 458          do iii=1,5
 459           iatom(iii)=iphilis(iii,i)+ishiftca
 460          enddo
 461
 462 C     ANGLE VECTOR CALCULTION
 463          RIX=C(1,IATOM(2))-C(1,IATOM(1))
 464          RIY=C(2,IATOM(2))-C(2,IATOM(1))
 465          RIZ=C(3,IATOM(2))-C(3,IATOM(1))
 466
 467          RIPX=C(1,IATOM(3))-C(1,IATOM(2))
 468          RIPY=C(2,IATOM(3))-C(2,IATOM(2))
 469          RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
 470
 471          RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
 472          RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
 473          RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
 474
 475          RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
 476          RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
 477          RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
 478
 479          GIX=RIY*RIPZ-RIZ*RIPY
 480          GIY=RIZ*RIPX-RIX*RIPZ
 481          GIZ=RIX*RIPY-RIY*RIPX
 482
 483          GIPX=RIPY*RIPPZ-RIPZ*RIPPY
 484          GIPY=RIPZ*RIPPX-RIPX*RIPPZ
 485          GIPZ=RIPX*RIPPY-RIPY*RIPPX
 486
 487          CIPX=C(1,IATOM(3))-C(1,IATOM(1))
 488          CIPY=C(2,IATOM(3))-C(2,IATOM(1))
 489          CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
 490
 491          CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
 492          CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
 493          CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
 494
 495          CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
 496          CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
 497          CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
 498
 499          DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
 500          DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
 501          DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
 502          DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
 503
 504 C     END OF ANGLE VECTOR CALCULTION
 505
 506          TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
 507          APHI(1)=TDOT/(DGI*DRIPP)
 508          TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
 509          APHI(2)=TDOT/(DGIP*DRIP3)
 510
 511          ephi = 0.0d0
 512          tfphi1=0.0d0
 513          tfphi2=0.0d0
 514          scc=0.0d0
 515
 516          do j=1, iphinum(i)
 517             DDPS1=APHI(1)-FPHI1(i,j)
 518             DDPS2=APHI(2)-FPHI2(i,j)
 519
 520             DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
 521
 522             if (dtmp.ge.15.0d0) then
 523                ps_tmp = 0.0d0
 524             else
 525 c               ps_tmp = dfaexp(idint(dtmp*1000)+1)
 526                 ps_tmp = exp(-dtmp)
 527             endif
 528
 529             ephi=ephi+sccphi(i,j)*ps_tmp
 530
 531             tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
 532             tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
 533
 534             scc=scc+sccphi(i,j)
 535 C            write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
 536 C     &           aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
 537          ENDDO
 538
 539          ephi=-ephi/scc*phi_inc*wwangle
 540          tfphi1=tfphi1/scc*phi_inc*wwangle
 541          tfphi2=tfphi2/scc*phi_inc*wwangle
 542
 543          IF (ABS(EPHI).LT.1d-20) THEN
 544             EPHI=0.0D0
 545          ENDIF
 546          IF (ABS(TFPHI1).LT.1d-20) THEN
 547             TFPHI1=0.0D0
 548          ENDIF
 549          IF (ABS(TFPHI2).LT.1d-20) THEN
 550             TFPHI2=0.0D0
 551          ENDIF
 552
 553 C     FORCE DIRECTION CALCULATION
 554          TDX(1:5)=0.0D0
 555          TDY(1:5)=0.0D0
 556          TDZ(1:5)=0.0D0
 557
 558          DM1=1.0d0/(DGI*DRIPP)
 559
 560          GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
 561          DM2=GIRPP/(DGI**3*DRIPP)
 562          DM3=GIRPP/(DGI*DRIPP**3)
 563
 564          DM4=1.0d0/(DGIP*DRIP3)
 565
 566          GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
 567          DM5=GIRP3/(DGIP**3*DRIP3)
 568          DM6=GIRP3/(DGIP*DRIP3**3)
 569 C     FIRST ATOM BY PHI1
 570          TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
 571      &        +( GIZ* RIPY- GIY* RIPZ)*DM2
 572          TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
 573      &        +( GIX* RIPZ- GIZ* RIPX)*DM2
 574          TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
 575      &        +( GIY* RIPX- GIX* RIPY)*DM2
 576          TDX(1)=TDX(1)*TFPHI1
 577          TDY(1)=TDY(1)*TFPHI1
 578          TDZ(1)=TDZ(1)*TFPHI1
 579 C     SECOND ATOM BY PHI1
 580          TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
 581      &        -(CIPY*GIZ-CIPZ*GIY)*DM2
 582          TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
 583      &        -(CIPZ*GIX-CIPX*GIZ)*DM2
 584          TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
 585      &        -(CIPX*GIY-CIPY*GIX)*DM2
 586          TDX(2)=TDX(2)*TFPHI1
 587          TDY(2)=TDY(2)*TFPHI1
 588          TDZ(2)=TDZ(2)*TFPHI1
 589 C     SECOND ATOM BY PHI2
 590          TDX(2)=TDX(2)+
 591      &        ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
 592      &        +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
 593          TDY(2)=TDY(2)+
 594      &        ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
 595      &        +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
 596          TDZ(2)=TDZ(2)+
 597      &        ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
 598      &        +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
 599 C     THIRD ATOM BY PHI1
 600          TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
 601      &        -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
 602          TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
 603      &        -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
 604          TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
 605      &        -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
 606          TDX(3)=TDX(3)*TFPHI1
 607          TDY(3)=TDY(3)*TFPHI1
 608          TDZ(3)=TDZ(3)*TFPHI1
 609 C     THIRD ATOM BY PHI2
 610          TDX(3)=TDX(3)+
 611      &        ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
 612      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
 613          TDY(3)=TDY(3)+
 614      &        ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
 615      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
 616          TDZ(3)=TDZ(3)+
 617      &        ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
 618      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
 619 C     FOURTH ATOM BY PHI1
 620          TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
 621          TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
 622          TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
 623 C     FOURTH ATOM BY PHI2
 624          TDX(4)=TDX(4)+
 625      &        ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
 626      &        -( GIPY*RIPZ-RIPY*GIPZ)*DM5
 627      &        + RIP3X*DM6)*TFPHI2
 628          TDY(4)=TDY(4)+
 629      &        ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
 630      &        -( GIPZ*RIPX-RIPZ*GIPX)*DM5
 631      &        + RIP3Y*DM6)*TFPHI2
 632          TDZ(4)=TDZ(4)+
 633      &        ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
 634      &        -( GIPX*RIPY-RIPX*GIPY)*DM5
 635      &        + RIP3Z*DM6)*TFPHI2
 636 C     FIFTH ATOM BY PHI2
 637          TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
 638          TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
 639          TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
 640 C     END OF FORCE DIRECTION
 641 c     force calcuation
 642          DO II=1,5
 643             gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
 644             gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
 645             gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
 646          ENDDO
 647 c     energy calculation
 648          enephi = enephi + ephi
 649 c     end of single assignment statement
 650       ENDDO
 651 C     END OF PHI RESTRAINT
 652
 653 C     START OF THETA ANGLE
 654       do i=idfathe_start,idfathe_end
 655
 656          athe = 0.0d0
 657          do iii=1,5
 658           iatom(iii)=ithelis(iii,i)+ishiftca
 659          enddo
 660
 661
 662 C     ANGLE VECTOR CALCULTION
 663          RIX=C(1,IATOM(2))-C(1,IATOM(1))
 664          RIY=C(2,IATOM(2))-C(2,IATOM(1))
 665          RIZ=C(3,IATOM(2))-C(3,IATOM(1))
 666
 667          RIPX=C(1,IATOM(3))-C(1,IATOM(2))
 668          RIPY=C(2,IATOM(3))-C(2,IATOM(2))
 669          RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
 670
 671          RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
 672          RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
 673          RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
 674
 675          RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
 676          RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
 677          RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
 678
 679          GIX=RIY*RIPZ-RIZ*RIPY
 680          GIY=RIZ*RIPX-RIX*RIPZ
 681          GIZ=RIX*RIPY-RIY*RIPX
 682
 683          GIPX=RIPY*RIPPZ-RIPZ*RIPPY
 684          GIPY=RIPZ*RIPPX-RIPX*RIPPZ
 685          GIPZ=RIPX*RIPPY-RIPY*RIPPX
 686
 687          GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
 688          GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
 689          GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
 690
 691          CIPX=C(1,IATOM(3))-C(1,IATOM(1))
 692          CIPY=C(2,IATOM(3))-C(2,IATOM(1))
 693          CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
 694
 695          CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
 696          CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
 697          CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
 698
 699          CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
 700          CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
 701          CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
 702
 703          DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
 704          DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
 705          DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
 706          DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
 707          DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
 708 C     END OF ANGLE VECTOR CALCULTION
 709
 710          TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
 711          ATHE(1)=TDOT/(DGI*DGIP)
 712          TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
 713          ATHE(2)=TDOT/(DGIP*DGIPP)
 714
 715          ETHE=0.0D0
 716          TFTHE1=0.0D0
 717          TFTHE2=0.0D0
 718          SCC=0.0D0
 719          TH_TMP=0.0d0
 720
 721          do j=1,ithenum(i)
 722             ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
 723             ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
 724             dtmp= (ddth1**2+ddth2**2)/cwidth2
 725             if ( dtmp .ge. 15.0d0) then
 726                th_tmp = 0.0d0
 727             else
 728 c               th_tmp = dfaexp ( idint(dtmp*1000)+1 )
 729                th_tmp = exp(-dtmp)
 730             end if
 731
 732             ethe=ethe+sccthe(i,j)*th_tmp
 733
 734             tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
 735             tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
 736             scc=scc+sccthe(i,j)
 737 c            write(2,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
 738 c     &           athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
 739          enddo
 740
 741          ethe=-ethe/scc*the_inc*wwangle
 742          tfthe1=tfthe1/scc*the_inc*wwangle
 743          tfthe2=tfthe2/scc*the_inc*wwangle
 744          IF (ABS(ETHE).LT.TENM20) THEN
 745             ETHE=0.0D0
 746          ENDIF
 747          IF (ABS(TFTHE1).LT.TENM20) THEN
 748             TFTHE1=0.0D0
 749          ENDIF
 750          IF (ABS(TFTHE2).LT.TENM20) THEN
 751             TFTHE2=0.0D0
 752          ENDIF
 753
 754          TDX(1:5)=0.0D0
 755          TDY(1:5)=0.0D0
 756          TDZ(1:5)=0.0D0
 757
 758          DM1=1.0d0/(DGI*DGIP)
 759          DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
 760          DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
 761
 762          DM4=1.0d0/(DGIP*DGIPP)
 763          DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
 764          DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
 765
 766 C     FIRST ATOM BY THETA1
 767          TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
 768      &        -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
 769          TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
 770      &        -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
 771          TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
 772      &        -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
 773 C     SECOND ATOM BY THETA1
 774          TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
 775      &        -(CIPY*GIZ-CIPZ*GIY)*DM2
 776      &        +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
 777          TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
 778      &        -(CIPZ*GIX-CIPX*GIZ)*DM2
 779      &        +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
 780          TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
 781      &        -(CIPX*GIY-CIPY*GIX)*DM2
 782      &        +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
 783 C     SECOND ATOM BY THETA2
 784          TDX(2)=TDX(2)+
 785      &        ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
 786      &        -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
 787          TDY(2)=TDY(2)+
 788      &        ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
 789      &        -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
 790          TDZ(2)=TDZ(2)+
 791      &        ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
 792      &        -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
 793 C     THIRD ATOM BY THETA1
 794          TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
 795      &        -(GIY*RIZ-GIZ*RIY)*DM2
 796      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
 797          TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
 798      &        -(GIZ*RIX-GIX*RIZ)*DM2
 799      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
 800          TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
 801      &        -(GIX*RIY-GIY*RIX)*DM2
 802      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
 803 C     THIRD ATOM BY THETA2
 804          TDX(3)=TDX(3)+
 805      &        ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
 806      &        -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
 807      &        +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
 808          TDY(3)=TDY(3)+
 809      &        ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
 810      &        -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
 811      &        +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
 812          TDZ(3)=TDZ(3)+
 813      &        ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
 814      &        -(CIPPX*GIPY-CIPPY*GIPX)*DM5
 815      &        +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
 816 C     FOURTH ATOM BY THETA1
 817          TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
 818      &        -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
 819          TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
 820      &        -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
 821          TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
 822      &        -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
 823 C     FOURTH ATOM BY THETA2
 824          TDX(4)=TDX(4)+
 825      &        ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
 826      &        -(GIPY*RIPZ-GIPZ*RIPY)*DM5
 827      &        -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
 828          TDY(4)=TDY(4)+
 829      &        ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
 830      &        -(GIPZ*RIPX-GIPX*RIPZ)*DM5
 831      &        -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
 832          TDZ(4)=TDZ(4)+
 833      &        ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
 834      &        -(GIPX*RIPY-GIPY*RIPX)*DM5
 835      &        -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
 836 C     FIFTH ATOM BY THETA2
 837          TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
 838      &        -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
 839          TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
 840      &        -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
 841          TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
 842      &        -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
 843 C     !! END OF FORCE DIRECTION!!!!
 844          DO II=1,5
 845             gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
 846             gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
 847             gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
 848          ENDDO
 849 C     energy calculation
 850          enethe = enethe + ethe
 851       ENDDO
 852
 853       edfator = enephi + enethe
 854
 855       RETURN
 856       END
 857
 858       subroutine edfan(edfanei)
 859 C     DFA neighboring CA restraint
 860       implicit real*8 (a-h,o-z)
 861       include 'DIMENSIONS'
 862       include 'COMMON.CHAIN'
 863       include 'COMMON.DERIV'
 864       include 'COMMON.DFA'
 865
 866       integer i,j,imin
 867       integer kshnum, n1atom
 868
 869       double precision enenei,tmp_n
 870       double precision pai,hpai
 871       double precision jix,jiy,jiz,ndiff,snorm_nei
 872       double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
 873       double precision dr,dr2,half,ntmp,dtmp
 874
 875       parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
 876       parameter(pai=3.14159265358979323846D0)
 877       parameter(hpai=1.5707963267948966D0)
 878       parameter(snorm_nei=0.886226925452758D0)
 879
 880       edfanei = 0.0d0
 881       enenei  = 0.0d0
 882       gdfan   = 0.0d0
 883
 884 c      print*, 's1:', s1(:)
 885 c      print*, 's2:', s2(:)
 886
 887       do i=idfanei_start,idfanei_end
 888
 889          kshnum=kshell(i)
 890          n1atom=ineilis(i)+ishiftca
 891 C         write(*,*) 'kshnum,n1atom:', kshnum, n1atom
 892
 893          tmp_n=0.0d0
 894          ftmp=0.0d0
 895          dnei=0.0d0
 896          dist=0.0d0
 897          t1dx=0.0d0
 898          t1dy=0.0d0
 899          t1dz=0.0d0
 900          t2dx=0.0d0
 901          t2dy=0.0d0
 902          t2dz=0.0d0
 903
 904          do j = ishiftca+1, ilastca
 905
 906             if (n1atom.eq.j) cycle
 907
 908             jix=c(1,j)-c(1,n1atom)
 909             jiy=c(2,j)-c(2,n1atom)
 910             jiz=c(3,j)-c(3,n1atom)
 911             dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
 912
 913 c            write(*,*) n1atom, j, dist
 914
 915             if(kshnum.ne.1)then
 916                if (dist.lt.s1(kshnum).and.
 917      &              dist.gt.s2(kshnum-1)) then
 918
 919                   tmp_n=tmp_n+1.0d0
 920
 921 c                  write(*,*) 'case1:',tmp_n
 922
 923 cc                  t1dx=t1dx+0.0d0
 924 cc                  t1dy=t1dy+0.0d0
 925 cc                  t1dz=t1dz+0.0d0
 926                   t2dx(j)=0.0d0
 927                   t2dy(j)=0.0d0
 928                   t2dz(j)=0.0d0
 929
 930                elseif(dist.ge.s1(kshnum).and.
 931      &                 dist.le.s2(kshnum)) then
 932
 933                   dnei=(dist-s1(kshnum))/dr2*pai
 934                   tmp_n=tmp_n + half*(1+cos(dnei))
 935 c                  write(*,*) 'case2:',tmp_n
 936                   ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
 937 c     center atom
 938                   t1dx=t1dx+jix*ftmp
 939                   t1dy=t1dy+jiy*ftmp
 940                   t1dz=t1dz+jiz*ftmp
 941 c     neighbor atoms
 942                   t2dx(j)=-jix*ftmp
 943                   t2dy(j)=-jiy*ftmp
 944                   t2dz(j)=-jiz*ftmp
 945 c
 946                elseif(dist.ge.s1(kshnum-1).and.
 947      &                 dist.le.s2(kshnum-1)) then
 948                   dnei=(dist-s1(kshnum-1))/dr2*pai
 949                   tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
 950 c                  write(*,*) 'case3:',tmp_n
 951                   ftmp = hpai*sin(dnei)/dr2/dist
 952 c     center atom
 953                   t1dx=t1dx+jix*ftmp
 954                   t1dy=t1dy+jiy*ftmp
 955                   t1dz=t1dz+jiz*ftmp
 956 c     neighbor atoms
 957                   t2dx(j)=-jix*ftmp
 958                   t2dy(j)=-jiy*ftmp
 959                   t2dz(j)=-jiz*ftmp
 960
 961                endif
 962
 963             elseif(kshnum.eq.1) then
 964
 965                if(dist.lt.s1(kshnum))then
 966
 967                   tmp_n=tmp_n+1.0d0
 968 c                  write(*,*) 'case4:',tmp_n
 969 cc                  t1dx=t1dx+0.0d0
 970 cc                  t1dy=t1dy+0.0d0
 971 cc                  t1dz=t1dz+0.0d0
 972                   t2dx(j)=0.0d0
 973                   t2dy(j)=0.0d0
 974                   t2dz(j)=0.0d0
 975
 976                elseif(dist.ge.s1(kshnum).and.
 977      &                 dist.le.s2(kshnum))then
 978
 979                   dnei=(dist-s1(kshnum))/dr2*pai
 980                   tmp_n=tmp_n + half*(1+cos(dnei))
 981 c                  write(*,*) 'case5:',tmp_n
 982                   ftmp = -hpai*sin(dnei)/dr2/dist
 983 c     center atom
 984                   t1dx=t1dx+jix*ftmp
 985                   t1dy=t1dy+jiy*ftmp
 986                   t1dz=t1dz+jiz*ftmp
 987 c     neighbor atoms
 988                   t2dx(j)=-jix*ftmp
 989                   t2dy(j)=-jiy*ftmp
 990                   t2dz(j)=-jiz*ftmp
 991
 992                endif
 993             endif
 994          enddo
 995
 996          scc=0.0d0
 997          enei=0.0d0
 998          tmp_fnei=0.0d0
 999          ndiff=0.0d0
1000
1001          do imin=1,ineinum(i)
1002
1003             ndiff = tmp_n-fnei(i,imin)
1004             dtmp  = ndiff*ndiff
1005
1006             if (dtmp.ge.15.0d0) then
1007                ntmp = 0.0d0
1008             else
1009 c               ntmp = dfaexp( idint(dtmp*1000) + 1 )
1010                 ntmp = exp(-dtmp)
1011             end if
1012
1013             enei=enei+sccnei(i,imin)*ntmp
1014             tmp_fnei=tmp_fnei-
1015      &           sccnei(i,imin)*ntmp*ndiff*2.0d0
1016             scc=scc+sccnei(i,imin)
1017
1018 c            write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
1019 c     &           fnei(i,imin),sccnei(i,imin),enei,scc
1020          enddo
1021
1022          enei=-enei/scc*snorm_nei*nei_inc*wwnei
1023          tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
1024
1025 c         if (abs(enei).lt.1.0d-20)then
1026 c            enei=0.0d0
1027 c         endif
1028 c         if (abs(tmp_fnei).lt.1.0d-20) then
1029 c            tmp_fnei=0.0d0
1030 c         endif
1031
1032 c     force calculation
1033          t1dx=t1dx*tmp_fnei
1034          t1dy=t1dy*tmp_fnei
1035          t1dz=t1dz*tmp_fnei
1036
1037          do j=ishiftca+1,ilastca
1038             t2dx(j)=t2dx(j)*tmp_fnei
1039             t2dy(j)=t2dy(j)*tmp_fnei
1040             t2dz(j)=t2dz(j)*tmp_fnei
1041          enddo
1042
1043          gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
1044          gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
1045          gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
1046
1047          do j=ishiftca+1,ilastca
1048             gdfan(1,j)=gdfan(1,j)+t2dx(j)
1049             gdfan(2,j)=gdfan(2,j)+t2dy(j)
1050             gdfan(3,j)=gdfan(3,j)+t2dz(j)
1051          enddo
1052 c     energy calculation
1053
1054          enenei=enenei+enei
1055
1056       enddo
1057
1058       edfanei=enenei
1059
1060       return
1061       end
1062
1063       subroutine edfab(edfabeta)
1064
1065       implicit real*8 (a-h,o-z)
1066
1067       include 'DIMENSIONS'
1068       include 'COMMON.CHAIN'
1069       include 'COMMON.DERIV'
1070       include 'COMMON.DFA'
1071
1072       real*8 PAI
1073       parameter(PAI=3.14159265358979323846D0)
1074       parameter (maxca=800)
1075 C     sheet variables
1076       real*8 bx(maxres),by(maxres),bz(maxres)
1077       real*8 vbet(maxres,maxres)
1078       real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
1079       real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
1080       real*8 vbeta,vbetp,vbetm
1081       real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
1082      &     c00,s00,ulnex,dnex
1083       real*8 dp45,dm45,w_beta
1084
1085       real*8 cph(maxca),cth(maxca)
1086       real*8 atx(maxca),aty(maxca),atz(maxca)
1087       real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1088       real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1089       real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1090       real*8 sth(maxca)
1091       real*8 astx(maxca),asty(maxca),astz(maxca)
1092       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1093       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1094       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1095
1096       real*8 atxnum(maxca),atynum(maxca),atznum(maxca),
1097      & astxnum(maxca),astynum(maxca),astznum(maxca),
1098      & atmxnum(maxca),atmynum(maxca),atmznum(maxca),
1099      & astmxnum(maxca),astmynum(maxca),astmznum(maxca),
1100      & atmmxnum(maxca),atmmynum(maxca),atmmznum(maxca),
1101      & astmmxnum(maxca),astmmynum(maxca),astmmznum(maxca),
1102      & atm3xnum(maxca),atm3ynum(maxca),atm3znum(maxca),
1103      & astm3xnum(maxca),astm3ynum(maxca),astm3znum(maxca),
1104      & cth_orig(maxca),sth_orig(maxca)
1105
1106       common /sheca/     bx,by,bz
1107       common /shee/      vbeta,vbet,vbetp,vbetm
1108       common /shetf/     shetfx,shetfy,shetfz
1109       common /shef/      shefx, shefy, shefz
1110       common /sheparm/   dca,dlhb,ulhb,dshe,dldhb,uldhb,
1111      &                   c00,s00,ulnex,dnex
1112       common /sheconst/  dp45,dm45,w_beta
1113
1114       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1115      $     atmmz,atm3x,atm3y,atm3z
1116       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1117      $     astmmz,astm3x,astm3y,astm3z
1118
1119       common /coscos/   cph,cth
1120       common /sinsin/ sth
1121
1122 C     End of sheet variables
1123
1124       integer i,j
1125       double precision enebet
1126
1127       enebet=0.0d0
1128 c      bx=0.0d0;by=0.0d0;bz=0.0d0
1129 c      shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1130
1131       gdfab=0.0d0
1132
1133       do i=ishiftca+1,ilastca
1134          bx(i-ishiftca)=c(1,i)
1135          by(i-ishiftca)=c(2,i)
1136          bz(i-ishiftca)=c(3,i)
1137       enddo
1138
1139 c      do i=1,ilastca-ishiftca
1140 c         read(99,*) bx(i),by(i),bz(i)
1141 c      enddo
1142 c      close(99)
1143
1144       dca=0.25d0**2
1145       dshe=0.3d0**2
1146       ULHB=5.0D0
1147       ULDHB=5.0D0
1148       ULNEX=COS(60.0D0/180.0D0*PAI)
1149
1150       DLHB=1.0D0
1151       DLDHB=1.0D0
1152
1153       DNEX=0.3D0**2
1154
1155       C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1156       S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1157
1158       W_BETA=0.5D0
1159       DP45=W_BETA
1160       DM45=W_BETA
1161
1162 C     END OF INITIALIZATION
1163
1164       nca=ilastca-ishiftca
1165
1166       call angvectors(nca)
1167       call sheetforce(nca,wshet)
1168
1169 c     end of sheet energy and force
1170
1171       do j=1,nca
1172          shetfx(j)=shetfx(j)*beta_inc
1173          shetfy(j)=shetfy(j)*beta_inc
1174          shetfz(j)=shetfz(j)*beta_inc
1175 c         write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1176       enddo
1177
1178       vbeta=vbeta*beta_inc
1179       enebet=vbeta
1180       edfabeta=enebet
1181
1182       do j=1,nca
1183          gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1184          gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1185          gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1186       enddo
1187
1188 #ifdef DEBUG1
1189       do j=1,nca
1190         write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1191       enddo
1192
1193
1194       gdfab=0
1195       dinc=0.001
1196       do j=1,nca
1197         cth_orig(j)=cth(j)
1198         sth_orig(j)=sth(j)
1199       enddo
1200
1201       do j=1,nca
1202
1203        bx(j)=bx(j)+dinc
1204        call angvectors(nca)
1205        bx(j)=bx(j)-2*dinc
1206        call angvectors(nca)
1207        atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1208        astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1209        if (j.gt.1) then
1210        atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1211        astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1212        endif
1213        if (j.gt.2) then
1214        atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1215        astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1216        endif
1217        if (j.gt.3) then
1218        atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1219        astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1220        endif
1221        bx(j)=bx(j)+dinc
1222        by(j)=by(j)+dinc
1223        call angvectors(nca)
1224        by(j)=by(j)-2*dinc
1225        call angvectors(nca)
1226        by(j)=by(j)+dinc
1227        atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1228        astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1229        if (j.gt.1) then
1230        atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1231        astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1232        endif
1233        if (j.gt.2) then
1234        atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1235        astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1236        endif
1237        if (j.gt.3) then
1238        atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1239        astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1240        endif
1241
1242        bz(j)=bz(j)+dinc
1243        call angvectors(nca)
1244        bz(j)=bz(j)-2*dinc
1245        call angvectors(nca)
1246        bz(j)=bz(j)+dinc
1247
1248        atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1249        astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1250        if (j.gt.1) then
1251        atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1252        astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1253        endif
1254        if (j.gt.2) then
1255        atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1256        astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1257        endif
1258        if (j.gt.3) then
1259        atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1260        astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1261        endif
1262
1263       enddo
1264
1265       do i=1,nca
1266         write (*,'(2i5,a2,6f10.5)')
1267      &  i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1268      &          astxnum(i),astx(i),astxnum(i)/astx(i),
1269      &  i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1270      &          astynum(i),asty(i),astynum(i)/asty(i),
1271      &  i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1272      &          astznum(i),astz(i),astznum(i)/astz(i),
1273      &  i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1274      &          astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1275      &  i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1276      &          astmynum(i),astmy(i),astmynum(i)/astmy(i),
1277      &  i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1278      &          astmznum(i),astmz(i),astmznum(i)/astmz(i),
1279      &  i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1280      &          astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1281      &  i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1282      &          astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1283      &  i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1284      &          astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1285      &  i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1286      &          astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1287      &  i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1288      &          astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1289      &  i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1290      &          astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1291      &  i,0," ",cth_orig(i),sth_orig(i)
1292       enddo
1293
1294
1295       gdfab=0
1296       dinc=0.001
1297
1298       do j=1,nca
1299
1300        bx(j)=bx(j)+dinc
1301        call angvectors(nca)
1302        call sheetforce(nca,wshet)
1303        vbeta1=vbeta*beta_inc
1304        bx(j)=bx(j)-2*dinc
1305        call angvectors(nca)
1306        call sheetforce(nca,wshet)
1307        vbeta2=vbeta*beta_inc
1308        gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1309        bx(j)=bx(j)+dinc
1310
1311        by(j)=by(j)+dinc
1312        call angvectors(nca)
1313        call sheetforce(nca,wshet)
1314        vbeta1=vbeta*beta_inc
1315        by(j)=by(j)-2*dinc
1316        call angvectors(nca)
1317        call sheetforce(nca,wshet)
1318        vbeta2=vbeta*beta_inc
1319        gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1320        by(j)=by(j)+dinc
1321
1322        bz(j)=bz(j)+dinc
1323        call angvectors(nca)
1324        call sheetforce(nca,wshet)
1325        vbeta1=vbeta*beta_inc
1326        bz(j)=bz(j)-2*dinc
1327        call angvectors(nca)
1328        call sheetforce(nca,wshet)
1329        vbeta2=vbeta*beta_inc
1330        gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1331        bz(j)=bz(j)+dinc
1332
1333
1334       enddo
1335
1336
1337       call angvectors(nca)
1338       call sheetforce(nca,wshet)
1339       do j=1,nca
1340          shetfx(j)=shetfx(j)*beta_inc
1341          shetfy(j)=shetfy(j)*beta_inc
1342          shetfz(j)=shetfz(j)*beta_inc
1343       enddo
1344
1345
1346       write(*,*) 'xyz analytical and numerical gradient'
1347       do j=1,nca
1348         write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1349      &                   ,(-gdfab(i,j),i=1,3)
1350       enddo
1351
1352       do j=1,nca
1353         write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1354      &                                  shetfy(j)/gdfab(2,j),
1355      &                                  shetfz(j)/gdfab(3,j)
1356       enddo
1357
1358       stop
1359 #endif
1360
1361       return
1362       end
1363 C-------------------------------------------------------------------------------
1364       subroutine angvectors(nca)
1365 c      implicit real*4(a-h,o-z)
1366       implicit none
1367       integer nca
1368       integer maxca
1369       parameter(maxca=800)
1370       real*8   pai,zero
1371       parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1372
1373       real*8   bx(maxca),by(maxca),bz(maxca)
1374       real*8   dis(maxca,maxca)
1375       real*8   apx(maxca),apy(maxca),apz(maxca)
1376       real*8   apmx(maxca),apmy(maxca),apmz(maxca)
1377       real*8   apmmx(maxca),apmmy(maxca),apmmz(maxca)
1378       real*8   apm3x(maxca),apm3y(maxca),apm3z(maxca)
1379       real*8   atx(maxca),aty(maxca),atz(maxca)
1380       real*8   atmx(maxca),atmy(maxca),atmz(maxca)
1381       real*8   atmmx(maxca),atmmy(maxca),atmmz(maxca)
1382       real*8   atm3x(maxca),atm3y(maxca),atm3z(maxca)
1383       real*8   astx(maxca),asty(maxca),astz(maxca)
1384       real*8   astmx(maxca),astmy(maxca),astmz(maxca)
1385       real*8   astmmx(maxca),astmmy(maxca),astmmz(maxca)
1386       real*8   astm3x(maxca),astm3y(maxca),astm3z(maxca)
1387       real*8   sth(maxca)
1388       real*8   cph(maxca),cth(maxca)
1389       real*8   ulcos(maxca)
1390       real*8   p,c
1391       integer  i, ip, ipp, ip3, j
1392       real*8   rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1393       real*8   rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1394       real*8   gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1395       real*8   cix, ciy, ciz, cipx, cipy, cipz
1396       real*8   gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1397       real*8   d10, d11, d12, d13, d20, d21, d22, d23, d24
1398       real*8   d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1399       real*8   d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1400       real*8   dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1401       real*8   dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1402       real*8   g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1403       real*8   gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1404       real*8   gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1405       real*8   gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1406       real*8   gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1407       real*8   grpp,gx,gy,gz
1408       real*8   rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1409       real*8   sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1410       integer inb,nmax,iselect
1411
1412       common /sheca/   bx,by,bz
1413       common /difvec/  rx, ry, rz
1414       common /ulang/    ulcos
1415       common /phys1/   inb,nmax,iselect
1416       common /phys4/   p,c
1417       common /kyori2/  dis
1418       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1419      &     apmmz,apm3x,apm3y,apm3z
1420       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1421      &     atmmz,atm3x,atm3y,atm3z
1422       common /coscos/   cph,cth
1423       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1424      &     astmmz,astm3x,astm3y,astm3z
1425       common /sinsin/   sth
1426 C-------------------------------------------------------------------------------
1427 c      write(*,*) 'inside angvectors'
1428 C     initialize
1429       p=0.1d0
1430       c=1.0d0
1431       inb=nca
1432       cph=zero; cth=zero; sth=zero
1433       apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1434       apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1435       atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1436       atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1437       astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1438       astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1439       astm3z=zero
1440 C     end of initialize
1441 C     r[x,y,z] calc and distance calculation
1442       rx=zero;ry=zero;rz=zero
1443
1444       do i=1,inb
1445          do j=1,inb
1446             rx(i,j)=bx(j)-bx(i)
1447             ry(i,j)=by(j)-by(i)
1448             rz(i,j)=bz(j)-bz(i)
1449             dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1450 c            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1451 c            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1452 c            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1453 c            write(*,*) 'dis(i,j):',i,j,dis(i,j)
1454          enddo
1455       enddo
1456 c     end of r[x,y,z] calc
1457 C     cos calc
1458       do i=1,inb-2
1459          ip=i+1
1460          ipp=i+2
1461
1462          if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1463             ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1464      $           +rz(i,ip)*rz(ip,ipp)
1465             ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1466          endif
1467       enddo
1468 c     end of virtual bond angle
1469 c      write(*,*) 'inside angvectors1'
1470 crc       do i=1,inb-3
1471       do i=1,inb
1472          ip=i+1
1473          ipp=i+2
1474          ip3=i+3
1475          rix=bx(ip)-bx(i)
1476          riy=by(ip)-by(i)
1477          riz=bz(ip)-bz(i)
1478          ripx=bx(ipp)-bx(ip)
1479          ripy=by(ipp)-by(ip)
1480          ripz=bz(ipp)-bz(ip)
1481          rippx=bx(ip3)-bx(ipp)
1482          rippy=by(ip3)-by(ipp)
1483          rippz=bz(ip3)-bz(ipp)
1484
1485          gx=riy*ripz-riz*ripy
1486          gy=riz*ripx-rix*ripz
1487          gz=rix*ripy-riy*ripx
1488          gpx=ripy*rippz-ripz*rippy
1489          gpy=ripz*rippx-ripx*rippz
1490          gpz=ripx*rippy-ripy*rippx
1491          gpcrp_x=gpy*ripz-gpz*ripy
1492          gpcrp_y=gpz*ripx-gpx*ripz
1493          gpcrp_z=gpx*ripy-gpy*ripx
1494          d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1495          gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1496      &        -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1497
1498          if(i.ge.2) then
1499             rimx=bx(i)-bx(i-1)
1500             rimy=by(i)-by(i-1)
1501             rimz=bz(i)-bz(i-1)
1502             gmx=rimy*riz-rimz*riy
1503             gmy=rimz*rix-rimx*riz
1504             gmz=rimx*riy-rimy*rix
1505             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1506             dgm3=dgm**3
1507             ggm=gmx*gx+gmy*gy+gmz*gz
1508             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1509             drim=dis(i-1,i)
1510             drim3=drim**3
1511             gcr_x=gy*riz-gz*riy
1512             gcr_y=gz*rix-gx*riz
1513             gcr_z=gx*riy-gy*rix
1514             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1515             d_gcr3=d_gcr**3
1516             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1517      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1518          endif
1519 c         write(*,*) 'inside angvectors2'
1520          if(i.ge.3) then
1521             rimmx=bx(i-1)-bx(i-2)
1522             rimmy=by(i-1)-by(i-2)
1523             rimmz=bz(i-1)-bz(i-2)
1524             drimm=dis(i-2,i-1)
1525             gmmx=rimmy*rimz-rimmz*rimy
1526             gmmy=rimmz*rimx-rimmx*rimz
1527             gmmz=rimmx*rimy-rimmy*rimx
1528             dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1529             dgmm3=dgmm**3
1530             gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1531             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1532             gmcrim_x=gmy*rimz-gmz*rimy
1533             gmcrim_y=gmz*rimx-gmx*rimz
1534             gmcrim_z=gmx*rimy-gmy*rimx
1535             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1536             d_gmcrim3=d_gmcrim**3
1537             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1538      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1539          endif
1540
1541          if(i.ge.4) then
1542             rim3x=bx(i-2)-bx(i-3)
1543             rim3y=by(i-2)-by(i-3)
1544             rim3z=bz(i-2)-bz(i-3)
1545             g3x=rim3y*rimmz-rim3z*rimmy
1546             g3y=rim3z*rimmx-rim3x*rimmz
1547             g3z=rim3x*rimmy-rim3y*rimmx
1548             dg30=sqrt(g3x**2+g3y**2+g3z**2)
1549             g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1550             g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1551 cc**********************************************************************
1552             gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1553             gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1554             gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1555             d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1556             d_gmmcrimm3=d_gmmcrimm**3
1557             gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1558      &           -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1559          endif
1560
1561          dri=dis(i,i+1)
1562          drip=dis(i+1,i+2)
1563          dripp=dis(i+2,i+3)
1564          dri3=dri**3
1565          dg=sqrt(gx**2+gy**2+gz**2)
1566          dgp=sqrt(gpx**2+gpy**2+gpz**2)
1567          dg3=dg**3
1568
1569          ggp=gx*gpx+gy*gpy+gz*gpz
1570          grpp=gx*rippx+gy*rippy+gz*rippz
1571
1572          if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1573      &        .and.d_gpcrp.gt.0.0D0) then
1574             cph(i)=grpp/dg/dripp
1575             cth(i)=ggp/dg/dgp
1576             sth(i)=gpcrp__g/d_gpcrp/dg
1577          else
1578 c
1579             cph(i)=1.0D0
1580             cth(i)=1.0D0
1581             sth(i)=0.0D0
1582          endif
1583
1584 c         write(*,*) 'inside angvectors3'
1585
1586          if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1587      &        .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1588             d10=1.0D0/(dg*dgp)
1589             d11=ggp/(dg3*dgp)
1590             d12=1.0D0/(dg*dripp)
1591             d13=grpp/(dg3*dripp)
1592             sd10=1.0D0/(d_gpcrp*dg)
1593             sd11=gpcrp__g/(d_gpcrp*dg3)
1594          else
1595             d10=0.0D0
1596             d11=0.0D0
1597             d12=0.0D0
1598             d13=0.0D0
1599             sd10=0.0D0
1600             sd11=0.0D0
1601          endif
1602
1603          atx(i)=(ripz*gpy-ripy*gpz)*d10
1604      &        -(gy*ripz-gz*ripy)*d11
1605          aty(i)=(ripx*gpz-ripz*gpx)*d10
1606      &        -(gz*ripx-gx*ripz)*d11
1607          atz(i)=(ripy*gpx-ripx*gpy)*d10
1608      &        -(gx*ripy-gy*ripx)*d11
1609          astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1610      &        +ripy*gpy*ripx-gpx*ripz**2)
1611      &        -sd11*(gy*ripz-gz*ripy)
1612          asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1613      &        -gpy*ripx**2+gpz*ripy*ripz)
1614      &        -sd11*(-gx*ripz+gz*ripx)
1615          astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1616      &        -gpz*ripy**2+ripz*gpx*ripx)
1617      &        -sd11*(gx*ripy-gy*ripx)
1618          apx(i)=(ripz*rippy-ripy*rippz)*d12
1619      &        -(gy*ripz-gz*ripy)*d13
1620          apy(i)=(ripx*rippz-ripz*rippx)*d12
1621      &        -(gz*ripx-gx*ripz)*d13
1622          apz(i)=(ripy*rippx-ripx*rippy)*d12
1623      &        -(gx*ripy-gy*ripx)*d13
1624
1625          if(i.ge.2) then
1626             cix=bx(ip)-bx(i-1)
1627             ciy=by(ip)-by(i-1)
1628             ciz=bz(ip)-bz(i-1)
1629             cipx=bx(ipp)-bx(i)
1630             cipy=by(ipp)-by(i)
1631             cipz=bz(ipp)-bz(i)
1632             ripx=bx(ipp)-bx(ip)
1633             ripy=by(ipp)-by(ip)
1634             ripz=bz(ipp)-bz(ip)
1635             if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1636      &           .and.d_gcr3.gt.0.0D0) then
1637                d20=1.0D0/(dg*dgm)
1638                d21=ggm/(dgm3*dg)
1639                d22=ggm/(dgm*dg3)
1640                d23=1.0D0/(dgm*drip)
1641                d24=gmrp/(dgm3*drip)
1642                sd20=1.0D0/(d_gcr*dgm)
1643                sd21=gcr__gm/(d_gcr3*dgm)
1644                sd22=gcr__gm/(d_gcr*dgm3)
1645             else
1646                d20=0.0D0
1647                d21=0.0D0
1648                d22=0.0D0
1649                d23=0.0D0
1650                d24=0.0D0
1651                sd20=0.0D0
1652                sd21=0.0D0
1653                sd22=0.0D0
1654             endif
1655             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1656      &           -(ciy*gmz-ciz*gmy)*d21
1657      &           +(ripy*gz-ripz*gy)*d22
1658             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1659      &           -(ciz*gmx-cix*gmz)*d21
1660      &           +(ripz*gx-ripx*gz)*d22
1661             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1662      &           -(cix*gmy-ciy*gmx)*d21
1663      &           +(ripx*gy-ripy*gx)*d22
1664 cc**********************************************************************
1665             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1666      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1667      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1668      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1669      &           +gcr_z*(-ripz*rix+gy))
1670      &           -sd22*(-gmy*ciz+gmz*ciy)
1671
1672             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1673      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1674      &           +riz*ripz*gmy)
1675      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1676      &           -gcr_z*(ripz*riy+gx))
1677      &           -sd22*(gmx*ciz-gmz*cix)
1678
1679             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1680      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1681      &           -riz*gx*cix)
1682      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1683      &           +gcr_z*(ripy*riy+ripx*rix))
1684      &           -sd22*(-gmx*ciy+gmy*cix)
1685 cc**********************************************************************
1686             apmx(i)=(ciy*ripz-ripy*ciz)*d23
1687      &           -(ciy*gmz-ciz*gmy)*d24
1688             apmy(i)=(ciz*ripx-ripz*cix)*d23
1689      &           -(ciz*gmx-cix*gmz)*d24
1690             apmz(i)=(cix*ripy-ripx*ciy)*d23
1691      &           -(cix*gmy-ciy*gmx)*d24
1692          endif
1693
1694          if(i.ge.3) then
1695             if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1696      &           .and.d_gmcrim3.gt.0.0D0) then
1697                d30=1.0D0/(dgm*dgmm)
1698                d31=gmmgm/(dgm3*dgmm)
1699                d32=gmmgm/(dgm*dgmm3)
1700                d33=1.0D0/(dgmm*dri)
1701                d34=gmmr/(dgmm3*dri)
1702                d35=gmmr/(dgmm*dri3)
1703                sd30=1.0D0/(d_gmcrim*dgmm)
1704                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1705                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1706             else
1707                d30=0.0D0
1708                d31=0.0D0
1709                d32=0.0D0
1710                d33=0.0D0
1711                d34=0.0D0
1712                d35=0.0D0
1713                sd30=0.0D0
1714                sd31=0.0D0
1715                sd32=0.0D0
1716             endif
1717
1718 c            write(*,*) 'inside angvectors4'
1719
1720 cc**********************************************************************
1721             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1722      &           -(ciy*gmz-ciz*gmy)*d31
1723      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1724             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1725      &           -(ciz*gmx-cix*gmz)*d31
1726      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1727             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1728      &           -(cix*gmy-ciy*gmx)*d31
1729      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1730 cc**********************************************************************
1731             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1732      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1733      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1734      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1735      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1736      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1737      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1738
1739             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1740      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1741      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1742      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1743      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1744      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1745      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1746
1747             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1748      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1749      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1750      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1751      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1752      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1753      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1754 c**********************************************************************
1755             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1756      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1757      &           +rix*d35
1758             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1759      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1760      &           +riy*d35
1761             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1762      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1763      &           +riz*d35
1764          endif
1765
1766          if(i.ge.4) then
1767             if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1768      &           .and.drim3.gt.0.0D0
1769      &           .and.d_gmmcrimm3.gt.0.0D0) then
1770                d40=1.0D0/(dg30*dgmm)
1771                d41=g3gmm/(dg30*dgmm3)
1772                d42=1.0D0/(dg30*drim)
1773                d43=g3rim_/(dg30*drim3)
1774                sd40=1.0D0/(dg30*d_gmmcrimm)
1775                sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1776             else
1777                d40=0.0D0
1778                d41=0.0D0
1779                d42=0.0D0
1780                d43=0.0D0
1781                sd40=0.0D0
1782                sd41=0.0D0
1783             endif
1784             atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1785      &           -(gmmy*rimmz-gmmz*rimmy)*d41
1786             atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1787      &           -(gmmz*rimmx-gmmx*rimmz)*d41
1788             atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1789      &           -(gmmx*rimmy-gmmy*rimmx)*d41
1790 cc**********************************************************************
1791             astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1792      &           -g3z*rimmz*rimmx+rimmy**2*g3x)
1793      &           -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1794      &           -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1795
1796             astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1797      &           -rimmx*rimmy*g3x+rimmz**2*g3y)
1798      &           -sd41*(-gmmcrimm_x*rimmx*rimmy
1799      &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1800
1801 c     &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1802
1803             astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1804      &           +g3z*rimmx**2-rimmz*rimmy*g3y)
1805      &           -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1806      &           +gmmcrimm_z*(rimmy**2+rimmx**2))
1807 c**********************************************************************
1808             apm3x(i)=g3x*d42-rimx*d43
1809             apm3y(i)=g3y*d42-rimy*d43
1810             apm3z(i)=g3z*d42-rimz*d43
1811          endif
1812       enddo
1813 c*******************************************************************************
1814
1815 c      write(*,*) 'inside angvectors5'
1816
1817 c       do i=inb-2,inb
1818        do i=1,0
1819          rimx=bx(i)-bx(i-1)
1820          rimy=by(i)-by(i-1)
1821          rimz=bz(i)-bz(i-1)
1822          rimmx=bx(i-1)-bx(i-2)
1823          rimmy=by(i-1)-by(i-2)
1824          rimmz=bz(i-1)-bz(i-2)
1825          rim3x=bx(i-2)-bx(i-3)
1826          rim3y=by(i-2)-by(i-3)
1827          rim3z=bz(i-2)-bz(i-3)
1828          gmmx=rimmy*rimz-rimmz*rimy
1829          gmmy=rimmz*rimx-rimmx*rimz
1830          gmmz=rimmx*rimy-rimmy*rimx
1831          g3x=rim3y*rimmz-rim3z*rimmy
1832          g3y=rim3z*rimmx-rim3x*rimmz
1833          g3z=rim3x*rimmy-rim3y*rimmx
1834
1835          dg30=sqrt(g3x**2+g3y**2+g3z**2)
1836          g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1837          dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1838          dgmm3=dgmm**3
1839          drim=dis(i-1,i)
1840          drimm=dis(i-2,i-1)
1841          drim3=drim**3
1842          g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1843 cc**********************************************************************
1844          gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1845          gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1846          gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1847          d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1848          d_gmmcrimm3=d_gmmcrimm**3
1849          gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1850      &        -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1851
1852          if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1853      &        .and.drim3.gt.0.0D0
1854      &        .and.d_gmmcrimm3.gt.0.0D0) then
1855             d40=1.0D0/(dg30*dgmm)
1856             d41=g3gmm/(dg30*dgmm3)
1857             d42=1.0D0/(dg30*drim)
1858             d43=g3rim_/(dg30*drim3)
1859             sd40=1.0D0/(dg30*d_gmmcrimm)
1860             sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1861          else
1862             d40=0.0D0
1863             d41=0.0D0
1864             d42=0.0D0
1865             d43=0.0D0
1866             sd40=0.0D0
1867             sd41=0.0D0
1868          endif
1869          atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1870      &        -(gmmy*rimmz-gmmz*rimmy)*d41
1871          atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1872      &        -(gmmz*rimmx-gmmx*rimmz)*d41
1873          atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1874      &        -(gmmx*rimmy-gmmy*rimmx)*d41
1875 cc**********************************************************************
1876          astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1877      &        -g3z*rimmz*rimmx+rimmy**2*g3x)
1878      &        -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1879      &        -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1880
1881          astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1882      &        -rimmx*rimmy*g3x+rimmz**2*g3y)
1883      &        -sd41*(-gmmcrimm_x*rimmx*rimmy
1884      &        +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1885
1886          astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1887      &        +g3z*rimmx**2-rimmz*rimmy*g3y)
1888      &        -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1889      &        +gmmcrimm_z*(rimmy**2+rimmx**2))
1890 cc**********************************************************************
1891          apm3x(i)=g3x*d42-rimx*d43
1892          apm3y(i)=g3y*d42-rimy*d43
1893          apm3z(i)=g3z*d42-rimz*d43
1894
1895          if(i.le.inb-1) then
1896             ip=i+1
1897             rix=bx(ip)-bx(i)
1898             riy=by(ip)-by(i)
1899             riz=bz(ip)-bz(i)
1900             cix=bx(ip)-bx(i-1)
1901             ciy=by(ip)-by(i-1)
1902             ciz=bz(ip)-bz(i-1)
1903             gmx=rimy*riz-rimz*riy
1904             gmy=rimz*rix-rimx*riz
1905             gmz=rimx*riy-rimy*rix
1906             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1907             dgm3=dgm**3
1908             dri=dis(i,i+1)
1909             dri3=dri**3
1910             gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1911             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1912             gmcrim_x=gmy*rimz-gmz*rimy
1913             gmcrim_y=gmz*rimx-gmx*rimz
1914             gmcrim_z=gmx*rimy-gmy*rimx
1915             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1916             d_gmcrim3=d_gmcrim**3
1917             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1918      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1919
1920             if(dgm3.gt.0.0D0.and.
1921      &           dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1922      &           .and.d_gmcrim3.gt.0.0D0) then
1923                d30=1.0D0/(dgm*dgmm)
1924                d31=gmmgm/(dgm3*dgmm)
1925                d32=gmmgm/(dgm*dgmm3)
1926                d33=1.0D0/(dgmm*dri)
1927                d34=gmmr/(dgmm3*dri)
1928                d35=gmmr/(dgmm*dri3)
1929                sd30=1.0D0/(d_gmcrim*dgmm)
1930                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1931                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1932
1933             else
1934                d30=0.0D0
1935                d31=0.0D0
1936                d32=0.0D0
1937                d33=0.0D0
1938                d34=0.0D0
1939                d35=0.0D0
1940                sd30=0.0D0
1941                sd31=0.0D0
1942                sd32=0.0D0
1943             endif
1944 cc**********************************************************************
1945             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1946      &           -(ciy*gmz-ciz*gmy)*d31
1947      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1948             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1949      &           -(ciz*gmx-cix*gmz)*d31
1950      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1951             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1952      &           -(cix*gmy-ciy*gmx)*d31
1953      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1954 cc**********************************************************************
1955             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1956      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1957      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1958      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1959      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1960      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1961      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1962
1963             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1964      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1965      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1966      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1967      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1968      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1969      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1970
1971             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1972      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1973      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1974      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1975      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1976      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1977      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1978 cc**********************************************************************
1979             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1980      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1981      &           +rix*d35
1982             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1983      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1984      &           +riy*d35
1985             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1986      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1987      &           +riz*d35
1988          endif
1989
1990 c         write(*,*) 'inside angvectors6'
1991
1992          if(i.eq.inb-2) then
1993             ipp=i+2
1994             ripx=bx(ipp)-bx(ip)
1995             ripy=by(ipp)-by(ip)
1996             ripz=bz(ipp)-bz(ip)
1997             cipx=bx(ipp)-bx(i)
1998             cipy=by(ipp)-by(i)
1999             cipz=bz(ipp)-bz(i)
2000             gx=riy*ripz-riz*ripy
2001             gy=riz*ripx-rix*ripz
2002             gz=rix*ripy-riy*ripx
2003             ggm=gmx*gx+gmy*gy+gmz*gz
2004             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
2005             dg=sqrt(gx**2+gy**2+gz**2)
2006             dg3=dg**3
2007             drip=dis(i+1,i+2)
2008             gcr_x=gy*riz-gz*riy
2009             gcr_y=gz*rix-gx*riz
2010             gcr_z=gx*riy-gy*rix
2011             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
2012             d_gcr3=d_gcr**3
2013             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
2014      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
2015             if(dgm3.gt.0.0D0.and.
2016      &           dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
2017      &           ) then
2018                d20=1.0D0/(dg*dgm)
2019                d21=ggm/(dgm3*dg)
2020                d22=ggm/(dgm*dg3)
2021                d23=1.0D0/(dgm*drip)
2022                d24=gmrp/(dgm3*drip)
2023                sd20=1.0D0/(d_gcr*dgm)
2024                sd21=gcr__gm/(d_gcr3*dgm)
2025                sd22=gcr__gm/(d_gcr*dgm3)
2026             else
2027                d20=0.0D0
2028                d21=0.0D0
2029                d22=0.0D0
2030                d23=0.0D0
2031                d24=0.0D0
2032                sd20=0.0D0
2033                sd21=0.0D0
2034                sd22=0.0D0
2035             endif
2036             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
2037      &           -(ciy*gmz-ciz*gmy)*d21
2038      &           +(ripy*gz-ripz*gy)*d22
2039             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
2040      &           -(ciz*gmx-cix*gmz)*d21
2041      &           +(ripz*gx-ripx*gz)*d22
2042             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
2043      &           -(cix*gmy-ciy*gmx)*d21
2044      &           +(ripx*gy-ripy*gx)*d22
2045 cc**********************************************************************
2046             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
2047      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
2048      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
2049      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
2050      &           +gcr_z*(-ripz*rix+gy))
2051      &           -sd22*(-gmy*ciz+gmz*ciy)
2052
2053             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
2054      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
2055      &           +riz*ripz*gmy)
2056      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
2057      &           -gcr_z*(ripz*riy+gx))
2058      &           -sd22*(gmx*ciz-gmz*cix)
2059
2060             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
2061      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
2062      &           -riz*gx*cix)
2063      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
2064      &           +gcr_z*(ripy*riy+ripx*rix))
2065      &           -sd22*(-gmx*ciy+gmy*cix)
2066 cc**********************************************************************
2067 c
2068             apmx(i)=(ciy*ripz-ripy*ciz)*d23
2069      &           -(ciy*gmz-ciz*gmy)*d24
2070             apmy(i)=(ciz*ripx-ripz*cix)*d23
2071      &           -(ciz*gmx-cix*gmz)*d24
2072             apmz(i)=(cix*ripy-ripx*ciy)*d23
2073      &           -(cix*gmy-ciy*gmx)*d24
2074
2075          endif
2076       enddo
2077
2078       return
2079       end
2080 c     END of angvectors
2081 c-------------------------------------------------------------------------------
2082 C---------------------------------------------------------------------------------
2083       subroutine sheetforce(nca,wshet)
2084       implicit none
2085 C     JYLEE
2086 c     this should be matched with dfa.fcm
2087       integer maxca
2088       parameter(maxca=800)
2089 cc**********************************************************************
2090       integer nca
2091       integer i,k
2092       integer inb,nmax,iselect
2093
2094 c      real*8 dfaexp(15001)
2095
2096       real*8 vbeta,vbetp,vbetm
2097       real*8 shefx(maxca,12)
2098       real*8 shefy(maxca,12),shefz(maxca,12)
2099       real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2100       real*8 vbet(maxca,maxca)
2101       real*8 wshet(maxca,maxca)
2102       real*8 bx(maxca),by(maxca),bz(maxca)
2103
2104       common /sheca/  bx,by,bz
2105       common /phys1/  inb,nmax,iselect
2106       common /shef/   shefx,shefy,shefz
2107       common /shee/   vbeta,vbet,vbetp,vbetm
2108       common /shetf/  shetfx,shetfy,shetfz
2109
2110       inb=nca
2111       do i=1,inb
2112          shetfx(i)=0.0D0
2113          shetfy(i)=0.0D0
2114          shetfz(i)=0.0D0
2115       enddo
2116
2117       do k=1,12
2118          do i=1,inb
2119             shefx(i,k)=0.0D0
2120             shefy(i,k)=0.0D0
2121             shefz(i,k)=0.0D0
2122          enddo
2123       enddo
2124
2125       call sheetene(nca,wshet)
2126       call sheetforce1
2127
2128  887  format(a,1x,i6,3x,f12.8)
2129  888  format(a,1x,i4,1x,i4,3x,f12.8)
2130  889  format(a,1x,i4,3x,f12.8)
2131       !write(2,*) 'coord : '
2132       do i=1,inb
2133          !write(2,887) 'bx:',i,bx(i)
2134          !write(2,887) 'by:',i,by(i)
2135          !write(2,887) 'bz:',i,bz(i)
2136       enddo
2137       !write(2,*) 'After sheetforce1'
2138       do i=1,inb
2139          do k=1,12
2140             !write(2,888) 'shefx :',i,k,shefx(i,k)
2141             !write(2,888) 'shefy :',i,k,shefy(i,k)
2142             !write(2,888) 'shefz :',i,k,shefz(i,k)
2143          enddo
2144       enddo
2145
2146       call sheetforce5
2147
2148       !write(2,*) 'After sheetforce5'
2149       do i=1,inb
2150          do k=1,12
2151             !write(2,888) 'shefx :',i,k,shefx(i,k)
2152             !write(2,888) 'shefy :',i,k,shefy(i,k)
2153             !write(2,888) 'shefz :',i,k,shefz(i,k)
2154          enddo
2155       enddo
2156
2157       call sheetforce6
2158
2159       !write(2,*) 'After sheetforce6'
2160       do i=1,inb
2161          do k=1,12
2162             !write(2,888) 'shefx :',i,k,shefx(i,k)
2163             !write(2,888) 'shefy :',i,k,shefy(i,k)
2164             !write(2,888) 'shefz :',i,k,shefz(i,k)
2165          enddo
2166       enddo
2167
2168       call sheetforce11
2169
2170       !write(2,*) 'After sheetforce11'
2171       do i=1,inb
2172          do k=1,12
2173             !write(2,888) 'shefx :',i,k,shefx(i,k)
2174             !write(2,888) 'shefy :',i,k,shefy(i,k)
2175             !write(2,888) 'shefz :',i,k,shefz(i,k)
2176          enddo
2177       enddo
2178
2179       call sheetforce12
2180
2181       !write(2,*) 'After sheetforce12'
2182       do i=1,inb
2183          do k=1,12
2184             !write(2,888) 'shefx :',i,k,shefx(i,k)
2185             !write(2,888) 'shefy :',i,k,shefy(i,k)
2186             !write(2,888) 'shefz :',i,k,shefz(i,k)
2187          enddo
2188       enddo
2189
2190       do i=1,inb
2191          do k=1,12
2192             shetfx(i)=shetfx(i)+shefx(i,k)
2193             shetfy(i)=shetfy(i)+shefy(i,k)
2194             shetfz(i)=shetfz(i)+shefz(i,k)
2195          enddo
2196       enddo
2197       !write(2,*) 'Beta Finished'
2198       do i=1,inb
2199          !write(2,889) 'shetfx : ',i,shetfx(i)
2200          !write(2,889) 'shetfy : ',i,shetfy(i)
2201          !write(2,889) 'shetfz : ',i,shetfz(i)
2202       enddo
2203
2204       return
2205       end
2206 C     end sheetforce
2207 c-------------------------------------------------------------------------------
2208       subroutine sheetene(nca,wshet)
2209       implicit none
2210       integer maxca
2211       parameter(maxca=800)
2212       real*8 dfa_cutoff,dfa_cutoff_delta
2213       parameter(dfa_cutoff=15.5d0)
2214       parameter(dfa_cutoff_delta=0.5d0)
2215 cc******************************************************************************
2216
2217 c      real*8 dfaexp(15001)
2218       real*8 dtmp1, dtmp2, dtmp3
2219
2220       real*8 vbet(maxca,maxca)
2221       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2222       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2223       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2224       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2225       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2226       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2227       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2228       real*8 cph(maxca),cth(maxca)
2229       real*8 rx(maxca,maxca)
2230       real*8 ry(maxca,maxca),rz(maxca,maxca)
2231       real*8 bx(maxca),by(maxca),bz(maxca)
2232       real*8 dis(maxca,maxca)
2233       real*8 ulcos(maxca)
2234 cc**********************************************************************
2235       real*8 astx(maxca),asty(maxca),astz(maxca)
2236       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2237       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2238       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2239       real*8 sth(maxca)
2240       real*8 wshet(maxca,maxca)
2241       real*8 dp45, dm45, w_beta
2242       real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2243       integer nca
2244       integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2245       real*8 uum, uup
2246       real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2247
2248       real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2249       common /shetf/  shetfx,shetfy,shetfz
2250
2251       common /sheca/    bx,by,bz
2252       common /phys1/    inb,nmax,iselect
2253       common /kyori2/   dis
2254       common /difvec/   rx,ry,rz
2255       common /coscos/   cph,cth
2256       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2257      &     c00,s00,ulnex,dnex
2258       common /sheconst/ dp45,dm45,w_beta
2259       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2260       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2261       common /shee/    vbeta,vbet,vbetp,vbetm
2262       common /ulang/    ulcos
2263 cc**********************************************************************
2264       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2265      &     astmmz,astm3x,astm3y,astm3z
2266       common /sinsin/   sth
2267
2268       real*8 r_pair_mat(maxca,maxca)
2269       real*8 e_gcont,fprim_gcont,de_gcont
2270 ci      integer istrand(maxca,maxca)
2271 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2272 ci      common  /shetest/ istrand,istrand_p,istrand_m
2273       common /beta_p/ r_pair_mat
2274 C-------------------------------------------------------------------------------
2275       r_pair_mat = 0.0d0
2276       do i=1,inb
2277          do j=1,inb
2278             r_pair_mat(i,j)=wshet(i,j)
2279 c            write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2280          enddo
2281       enddo
2282 c      stop
2283 c
2284       vbeta=0.0D0
2285       vbetp=0.0D0
2286       vbetm=0.0D0
2287
2288       do i=1,inb-7
2289          do j=i+4,inb-3
2290
2291             if (dis(i,j).lt.dfa_cutoff) then
2292             call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2293      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
2294
2295
2296             ip=i+1
2297             ipp=i+2
2298             jp=j+1
2299             jpp=j+2
2300 cc**********************************************************************
2301             y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2302      &           +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2303             y1=-0.5d0*y1/dca
2304             y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2305      &           +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2306             y2=-0.5d0*y2/dnex
2307
2308 cdebug            y2=0
2309
2310             y=y1+y2
2311
2312 ci           if(y.ge.-4) then
2313 ci              istrand(i,j)=1
2314 ci           else
2315 ci              istrand(i,j)=0
2316 ci           endif
2317
2318 ci           if(istrand(i,j).eq.1) then
2319
2320             yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2321             yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2322
2323
2324             pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2325      $           +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2326             pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2327      $           +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2328             pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2329      $           +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2330             pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2331      $           +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2332
2333            yshe1=pin1(i,j)**2+pin2(i,j)**2
2334            yshe1=-0.5d0*yshe1/dshe
2335            yshe2=pin3(i,j)**2+pin4(i,j)**2
2336            yshe2=-0.5d0*yshe2/dshe
2337
2338 ci              if((yshe1+yshe2).ge.-4) then
2339 ci                 istrand_p(i,j)=1
2340 ci              else
2341 ci                 istrand_p(i,j)=0
2342 ci              endif
2343
2344
2345 C            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2346 C            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2347 C            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2348 C            write(*,*) 'dis(i,j):',i,j,dis(i,j)
2349 C            write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2350 C            write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2351 C            write(*,*) 'pin1:',pin1(i,j)
2352 C            write(*,*) 'pin2:',pin2(i,j)
2353 C            write(*,*) 'pin3:',pin3(i,j)
2354 C            write(*,*) 'pin4:',pin4(i,j)
2355
2356 C            write(*,*) 'y:',y
2357 C            write(*,*) 'yy1:',yy1
2358 C            write(*,*) 'yy2:',yy2
2359 C            write(*,*) 'yshe1:',yshe1
2360 C            write(*,*) 'yshe2:',yshe2
2361 c
2362
2363 ci           if (istrand_p(i,j).eq.1) then
2364
2365 cd           yy1=0
2366 cd           yy2=0
2367 cd           yshe1=0
2368 cd           yshe2=0
2369            dtmp1 = y+yy1+yshe1
2370            dtmp2 = y+yy2+yshe2
2371            dtmp3 = y+yy1+yy2+yshe1+yshe2
2372
2373 C            write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2374 C            write(*,*)'2', y,yy1,yy2
2375 C            write(*,*)'3', yshe1,yshe2
2376
2377 cc           if (dtmp3.le.-35.0d0) then
2378 c              vbetap(i,j)=-dp45*exp(dtmp3)
2379 cc              vbetap(i,j)=0.0d0
2380 cc           else
2381 c              vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2382               vbetap(i,j)=-dp45*exp(dtmp3)
2383 cc           end if
2384
2385 cc           if (dtmp1.le.-35.0d0) then
2386 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2387 cc              vbetap1(i,j)=0.0d0
2388 cc           else
2389 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2390 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2391                vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2392 cc           end if
2393
2394 cc           if (dtmp2.le.-35.0d0) then
2395 C              vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2396 cc              vbetap2(i,j)=0.0d0
2397 cc           else
2398 c              vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2399 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2400               vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2401 cc           end if
2402
2403 c           vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2404 c           vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2405 c           vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2406
2407 !           write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2408 !           write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2409
2410 ci           elseif (istrand_p(i,j).eq.0)then
2411 ci            vbetap(i,j)=0
2412 ci            vbetap1(i,j)=0
2413 ci            vbetap2(i,j)=0
2414 ci           endif
2415
2416            yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2417            yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2418
2419            pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2420      $          +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2421            pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2422      $          +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2423            pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2424      $          +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2425            pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2426      $          +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2427
2428            yshe1=pina1(i,j)**2+pina2(i,j)**2
2429            yshe1=-0.5d0*yshe1/dshe
2430            yshe2=pina3(i,j)**2+pina4(i,j)**2
2431            yshe2=-0.5d0*yshe2/dshe
2432
2433 ci              if((yshe1+yshe2).ge.-4) then
2434 ci                 istrand_m(i,j)=1
2435 ci              else
2436 ci                 istrand_m(i,j)=0
2437 ci              endif
2438
2439
2440 C            write(*,*) 'pina1:',pina1(i,j)
2441 C            write(*,*) 'pina2:',pina2(i,j)
2442 C            write(*,*) 'pina3:',pina3(i,j)
2443 C            write(*,*) 'pina4:',pina4(i,j)
2444 C            write(*,*) 'yshe1:',yshe1
2445 C            write(*,*) 'yshe2:',yshe2
2446 C            write(*,*) 'dshe:',dshe
2447
2448 ci           if (istrand_m(i,j).eq.1) then
2449
2450 cd           yy1=0
2451 cd           yy2=0
2452 cd           yshe1=0
2453 cd           yshe2=0
2454
2455            dtmp3=y+yy1+yy2+yshe1+yshe2
2456            dtmp1=y+yy1+yshe1
2457            dtmp2=y+yy2+yshe2
2458
2459 cc           if(dtmp3 .le. -35.0d0) then
2460 c              vbetam(i,j)=-dm45*exp(dtmp3)
2461 cc              vbetam(i,j)=0.0d0
2462 cc           else
2463 c              vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2464               vbetam(i,j)=-dm45*exp(dtmp3)
2465 cc           end if
2466
2467 cc           if(dtmp1 .le. -35.0d0) then
2468 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2469 cc               vbetam1(i,j)=0.0d0
2470 cc           else
2471 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2472 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2473                vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2474 cc           end if
2475
2476 cc           if(dtmp2.le.-35.0d0) then
2477 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2478 cc              vbetam2(i,j)=0.0d0
2479 cc           else
2480 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2481 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2482               vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2483 cc           end if
2484
2485 ci           elseif (istrand_m(i,j).eq.0)then
2486 ci            vbetam(i,j)=0
2487 ci            vbetam1(i,j)=0
2488 ci            vbetam2(i,j)=0
2489 ci           endif
2490
2491
2492 c           vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2493 c           vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2494 c           vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2495
2496 !           write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2497 !           write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2498
2499            uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2500            uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2501
2502 c           write(*,*) 'uup,uum:', uup, uum
2503
2504 c           uup=vbetap1(i,j)+vbetap2(i,j)
2505 c           uum=vbetam1(i,j)+vbetam2(i,j)
2506
2507            vbet(i,j)=uup+uum
2508            vbetp=vbetp+uup
2509            vbetm=vbetm+uum
2510            vbeta=vbeta+vbet(i,j)*e_gcont
2511
2512
2513            if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
2514 c gradient correction from gcont
2515              de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
2516              shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
2517              shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
2518              shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
2519
2520              shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
2521              shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
2522              shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
2523
2524 c energy correction from gcont
2525              vbet(i,j)=vbet(i,j)*e_gcont
2526              vbetap(i,j)=vbetap(i,j)*e_gcont
2527              vbetap1(i,j)=vbetap1(i,j)*e_gcont
2528              vbetap2(i,j)=vbetap2(i,j)*e_gcont
2529              vbetam(i,j)=vbetam(i,j)*e_gcont
2530              vbetam1(i,j)=vbetam1(i,j)*e_gcont
2531              vbetam2(i,j)=vbetam2(i,j)*e_gcont
2532            endif
2533
2534
2535 ci         elseif(istrand(i,j).eq.0)then
2536 ci           vbet(i,j)=0
2537 ci         endif
2538
2539 c           write(*,*) 'uup,uum:',uup,uum
2540 c           write(*,*) 'vbetap(i,j):',vbetap(i,j)
2541 c           write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2542 c           write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2543 c           write(*,*) 'vbetam(i,j):',vbetam(i,j)
2544 c           write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2545 c           write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2546 c           write(*,*) 'uup:',uup
2547 c           write(*,*) 'uum:',uum
2548 c           write(*,*) 'vbetp:',vbetp
2549 c           write(*,*) 'vbetm:',vbetm
2550 c           write(*,*) 'vbet(i,j):',vbet(i,j)
2551 c           stop
2552
2553             else
2554               vbetap(i,j)=0
2555               vbetap1(i,j)=0
2556               vbetap2(i,j)=0
2557               vbetam(i,j)=0
2558               vbetam1(i,j)=0
2559               vbetam2(i,j)=0
2560               vbet(i,j)=0
2561             endif
2562         enddo
2563       enddo
2564
2565 !      do i=1,inb-7
2566 !         do j=i+4,inb-3
2567 !            write(*,*) 'I,J:', i,j
2568 !            write(*,*) 'vbetap(i,j):',vbetap(i,j)
2569 !            write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2570 !            write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2571 !            write(*,*) 'vbetam(i,j):',vbetam(i,j)
2572 !            write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2573 !            write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2574 !            write(*,*) 'vbet(i,j):',vbet(i,j)
2575 !         enddo
2576 !      enddo
2577
2578       return
2579       end
2580 c-------------------------------------------------------------------------------
2581       subroutine sheetforce1
2582       implicit none
2583       integer maxca
2584       parameter(maxca=800)
2585       real*8 dfa_cutoff,dfa_cutoff_delta
2586       parameter(dfa_cutoff=15.5d0)
2587       parameter(dfa_cutoff_delta=0.5d0)
2588 cc**********************************************************************
2589       real*8 vbet(maxca,maxca)
2590       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2591       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2592       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2593       real*8 cph(maxca),cth(maxca)
2594       real*8 rx(maxca,maxca)
2595       real*8 ry(maxca,maxca),rz(maxca,maxca)
2596       real*8 bx(maxca),by(maxca),bz(maxca)
2597       real*8 dis(maxca,maxca)
2598       real*8 shefx(maxca,12)
2599       real*8 shefy(maxca,12),shefz(maxca,12)
2600       real*8 atx(maxca),aty(maxca),atz(maxca)
2601       real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2602       real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2603       real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2604       real*8 apx(maxca),apy(maxca),apz(maxca)
2605       real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2606       real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2607       real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2608       real*8 ulcos(maxca)
2609       real*8 astx(maxca),asty(maxca),astz(maxca)
2610       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2611       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2612       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2613       real*8 sth(maxca)
2614       real*8 w_beta,dp45, dm45
2615       real*8 vbeta, vbetp, vbetm
2616       real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2617      $     c00,s00,ulnex,dnex
2618       integer inb,nmax,iselect
2619
2620       common /phys1/     inb,nmax,iselect
2621       common /kyori2/    dis
2622       common /difvec/   rx,ry,rz
2623       common /coscos/   cph,cth
2624       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2625      $     c00,s00,ulnex,dnex
2626       common /sheconst/ dp45,dm45,w_beta
2627       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2628       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2629      $     atmmz,atm3x,atm3y,atm3z
2630       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2631      $     apmmz,apm3x,apm3y,apm3z
2632       common /shef/   shefx,shefy,shefz
2633       common /shee/   vbeta,vbet,vbetp,vbetm
2634       common /ulang/    ulcos
2635 c     c**********************************************************************
2636       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2637      $     astmmz,astm3x,astm3y,astm3z
2638       common /sinsin/   sth
2639 C--------------------------------------------------------------------------------
2640 c     local variables
2641       integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2642       real*8  c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2643       real*8  c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2644       real*8  c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2645       real*8  dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2646       real*8 e_gcont,fprim_gcont
2647 C--------------------------------------------------------------------------------
2648       do i=4,inb-4
2649          im3=i-3
2650          imm=i-2
2651          im=i-1
2652          c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2653          v1=0.0D0
2654          do j=i+1,inb-3
2655             v1=v1+vbet(im3,j)
2656          enddo
2657          cc1=(ulcos(imm)-ulnex)/dnex
2658          dmm=cc1/(dis(imm,im)*dis(im,i))
2659          dmm__=cc1*ulcos(imm)/dis(im,i)**2
2660          fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2661          fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2662          fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2663 cd         fx=0
2664 cd         fy=0
2665 cd         fz=0
2666          fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2667          fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2668          fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2669          shefx(i,1)=fx*v1
2670          shefy(i,1)=fy*v1
2671          shefz(i,1)=fz*v1
2672       enddo
2673
2674       do i=3,inb-5
2675          imm=i-2
2676          im=i-1
2677          ip=i+1
2678          c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2679          v2=0.0D0
2680          do j=i+2,inb-3
2681             v2=v2+vbet(imm,j)
2682          enddo
2683          cc1=(ulcos(imm)-ulnex)/dnex
2684          cc2=(ulcos(im)-ulnex)/dnex
2685          dmm1=cc1/(dis(imm,im)*dis(im,i))
2686          dmm2=cc2/(dis(im,i)*dis(i,ip))
2687          dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2688          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2689          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2690 cc**********************************************************************
2691          fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2692      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2693          fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2694      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2695          fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2696      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2697 cd         fx=0
2698 cd         fy=0
2699 cd         fz=0
2700          fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2701          fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2702          fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2703          shefx(i,2)=fx*v2
2704          shefy(i,2)=fy*v2
2705          shefz(i,2)=fz*v2
2706       enddo
2707       do i=2,inb-6
2708          im=i-1
2709          ip=i+1
2710          ipp=i+2
2711          c3=(cth(im)*c00+sth(im)*s00-1)/dca
2712          v3=0.0D0
2713          do j=i+3,inb-3
2714             v3=v3+vbet(im,j)
2715          enddo
2716          cc2=(ulcos(im)-ulnex)/dnex
2717          cc3=(ulcos(i)-ulnex)/dnex
2718          dmm2=cc2/(dis(im,i)*dis(i,ip))
2719          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2720          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2721          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2722          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2723          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2724      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2725          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2726      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2727          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2728      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2729 cd         fx=0
2730 cd         fy=0
2731 cd         fz=0
2732          fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2733          fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2734          fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2735          shefx(i,3)=fx*v3
2736          shefy(i,3)=fy*v3
2737          shefz(i,3)=fz*v3
2738       enddo
2739       do i=1,inb-7
2740          ip=i+1
2741          ipp=i+2
2742          c4=(cth(i)*c00+sth(i)*s00-1)/dca
2743          v4=0.0D0
2744          do j=i+4,inb-3
2745             v4=v4+vbet(i,j)
2746          enddo
2747          cc3=(ulcos(i)-ulnex)/dnex
2748          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2749          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2750          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2751          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2752          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2753 cd         fx=0
2754 cd         fy=0
2755 cd         fz=0
2756          fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2757          fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2758          fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2759          shefx(i,4)=fx*v4
2760          shefy(i,4)=fy*v4
2761          shefz(i,4)=fz*v4
2762       enddo
2763       do j=8,inb
2764          jm3=j-3
2765          jmm=j-2
2766          jm=j-1
2767          c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2768          v7=0.0D0
2769          do i=1,j-7
2770             v7=v7+vbet(i,jm3)
2771          enddo
2772          cc7=(ulcos(jmm)-ulnex)/dnex
2773          dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2774          dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2775          fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2776          fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2777          fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2778 cd         fx=0
2779 cd         fy=0
2780 cd         fz=0
2781          fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2782          fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2783          fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2784          shefx(j,7)=fx*v7
2785          shefy(j,7)=fy*v7
2786          shefz(j,7)=fz*v7
2787       enddo
2788       do j=7,inb-1
2789          jm=j-1
2790          jmm=j-2
2791          jp=j+1
2792          c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2793          v8=0.0D0
2794          do i=1,j-6
2795             v8=v8+vbet(i,jmm)
2796          enddo
2797          cc7=(ulcos(jmm)-ulnex)/dnex
2798          cc8=(ulcos(jm)-ulnex)/dnex
2799          dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2800          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2801          dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2802          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2803          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2804          fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2805      $        -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2806          fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2807      $        -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2808          fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2809      $        -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2810 cd         fx=0
2811 cd         fy=0
2812 cd         fz=0
2813          fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2814          fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2815          fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2816          shefx(j,8)=fx*v8
2817          shefy(j,8)=fy*v8
2818          shefz(j,8)=fz*v8
2819       enddo
2820
2821       do j=6,inb-2
2822          jm=j-1
2823          jp=j+1
2824          jpp=j+2
2825          c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2826          v9=0.0D0
2827          do i=1,j-5
2828             v9=v9+vbet(i,jm)
2829          enddo
2830          cc8=(ulcos(jm)-ulnex)/dnex
2831          cc9=(ulcos(j)-ulnex)/dnex
2832          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2833          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2834          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2835          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2836          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2837          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2838      $        -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2839          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2840      $        -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2841          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2842      $        -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2843 cd         fx=0
2844 cd         fy=0
2845 cd         fz=0
2846          fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2847          fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2848          fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2849          shefx(j,9)=fx*v9
2850          shefy(j,9)=fy*v9
2851          shefz(j,9)=fz*v9
2852       enddo
2853
2854       do j=5,inb-3
2855          jp=j+1
2856          jpp=j+2
2857          c10=(cth(j)*c00+sth(j)*s00-1)/dca
2858          v10=0.0D0
2859          do i=1,j-4
2860             v10=v10+vbet(i,j)
2861          enddo
2862          cc9=(ulcos(j)-ulnex)/dnex
2863          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2864          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2865          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2866          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2867          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2868 cd         fx=0
2869 cd         fy=0
2870 cd         fz=0
2871          fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2872          fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2873          fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2874          shefx(j,10)=fx*v10
2875          shefy(j,10)=fy*v10
2876          shefz(j,10)=fz*v10
2877       enddo
2878
2879       return
2880       end
2881 c----------------------------------------------------------------------------
2882       subroutine sheetforce5
2883       implicit none
2884       integer maxca
2885       parameter(maxca=800)
2886       real*8 dfa_cutoff,dfa_cutoff_delta
2887       parameter(dfa_cutoff=15.5d0)
2888       parameter(dfa_cutoff_delta=0.5d0)
2889 cc**********************************************************************
2890       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2891       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2892       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2893       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2894       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2895       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2896       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2897       real*8 rx(maxca,maxca)
2898       real*8 ry(maxca,maxca),rz(maxca,maxca)
2899       real*8 bx(maxca),by(maxca),bz(maxca)
2900       real*8 dis(maxca,maxca)
2901       real*8 shefx(maxca,12),shefy(maxca,12)
2902       real*8 shefz(maxca,12)
2903       real*8 dp45,dm45,w_beta
2904       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2905      $     c00,s00,ulnex,dnex
2906       integer    inb,nmax,iselect
2907 cc**********************************************************************
2908       common /phys1/     inb,nmax,iselect
2909       common /kyori2/    dis
2910       common /difvec/   rx,ry,rz
2911       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2912      $     c00,s00,ulnex,dnex
2913       common /sheconst/ dp45,dm45,w_beta
2914       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2915       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2916       common /shef/   shefx,shefy,shefz
2917 ci      integer istrand(maxca,maxca)
2918 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2919 ci      common  /shetest/ istrand,istrand_p,istrand_m
2920 c********************************************************************************
2921 c     local variables
2922       integer i,imm,im,jp,jpp,j
2923       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2924       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2925       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2926       real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2927       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2928       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2929       real*8 e_gcont,fprim_gcont
2930 c********************************************************************************
2931       do i=3,inb-5
2932          imm=i-2
2933          im=i-1
2934          do j=i+2,inb-3
2935
2936             if (dis(imm,j).lt.dfa_cutoff) then
2937             call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
2938      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
2939
2940             jp=j+1
2941             jpp=j+2
2942
2943 ci            if(istrand(imm,j).eq.1
2944 ci     &   .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2945
2946
2947             yy1=-(dis(i,jpp)-ulhb)/dlhb
2948             y1x=rx(jpp,i)/dis(i,jpp)
2949             y1y=ry(jpp,i)/dis(i,jpp)
2950             y1z=rz(jpp,i)/dis(i,jpp)
2951             y11x=yy1*y1x
2952             y11y=yy1*y1y
2953             y11z=yy1*y1z
2954
2955             yy33=1.0D0/(dis(im,jp)*dis(im,i))
2956             yyy3=pin1(imm,j)/(dis(im,i)**2)
2957             yy3=-pin1(imm,j)/dshe
2958             y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2959             y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2960             y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2961
2962             yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2963             yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2964             yyy4b=pin3(imm,j)/(dis(im,i)**2)
2965             yy4=-pin3(imm,j)/dshe
2966             y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2967      $           -yyy4b*rx(im,i))*yy4
2968             y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2969      $           -yyy4b*ry(im,i))*yy4
2970             y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2971      $           -yyy4b*rz(im,i))*yy4
2972
2973
2974             yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2975             yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2976             yy5=-pin4(imm,j)/dshe
2977             y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2978             y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2979             y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2980
2981             sx=y11x+y3x+y4x+y5x
2982             sy=y11y+y3y+y4y+y5y
2983             sz=y11z+y3z+y4z+y5z
2984
2985             sx1=y3x
2986             sy1=y3y
2987             sz1=y3z
2988             sx2=y11x+y4x+y5x
2989             sy2=y11y+y4y+y5y
2990             sz2=y11z+y4z+y5z
2991
2992             shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2993      $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2994             shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2995      $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2996             shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2997      $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2998
2999 !            shefx(i,5)=shefx(i,5)
3000 !     $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
3001 !            shefy(i,5)=shefy(i,5)
3002 !     $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
3003 !            shefz(i,5)=shefz(i,5)
3004 !     $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
3005
3006             yy6=-(dis(i,jp)-uldhb)/dldhb
3007             y6x=rx(jp,i)/dis(i,jp)
3008             y6y=ry(jp,i)/dis(i,jp)
3009             y6z=rz(jp,i)/dis(i,jp)
3010             y66x=yy6*y6x
3011             y66y=yy6*y6y
3012             y66z=yy6*y6z
3013
3014             yy88=1.0D0/(dis(im,jpp)*dis(im,i))
3015             yyy8=pina1(imm,j)/(dis(im,i)**2)
3016             yy8=-pina1(imm,j)/dshe
3017             y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
3018             y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
3019             y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
3020
3021             yy99=1.0D0/(dis(jp,i)*dis(im,i))
3022             yyy9a=pina3(imm,j)/(dis(jp,i)**2)
3023             yyy9b=pina3(imm,j)/(dis(im,i)**2)
3024             yy9=-pina3(imm,j)/dshe
3025             y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
3026      $           -yyy9b*rx(im,i))*yy9
3027             y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
3028      $           -yyy9b*ry(im,i))*yy9
3029             y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
3030      $           -yyy9b*rz(im,i))*yy9
3031
3032             yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
3033             yyy10=pina4(imm,j)/(dis(jp,i)**2)
3034             yy10=-pina4(imm,j)/dshe
3035             y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
3036             y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
3037             y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
3038
3039             sx=y66x+y8x+y9x+y10x
3040             sy=y66y+y8y+y9y+y10y
3041             sz=y66z+y8z+y9z+y10z
3042
3043             sx1=y8x
3044             sy1=y8y
3045             sz1=y8z
3046             sx2=y66x+y9x+y10x
3047             sy2=y66y+y9y+y10y
3048             sz2=y66z+y9z+y10z
3049
3050             shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
3051      $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3052            shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
3053      $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3054             shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
3055      $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3056
3057 !            shefx(i,5)=shefx(i,5)
3058 !     $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
3059 !            shefy(i,5)=shefy(i,5)
3060 !     $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
3061 !            shefz(i,5)=shefz(i,5)
3062 !     $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
3063
3064              endif
3065 ci          endif
3066
3067          enddo
3068       enddo
3069
3070       return
3071       end
3072 c--------------------------------------------------------------------------c
3073       subroutine sheetforce6
3074       implicit none
3075       integer maxca
3076       parameter(maxca=800)
3077       real*8 dfa_cutoff,dfa_cutoff_delta
3078       parameter(dfa_cutoff=15.5d0)
3079       parameter(dfa_cutoff_delta=0.5d0)
3080 cc**********************************************************************
3081       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3082       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3083       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3084       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3085       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3086       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3087       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3088       real*8 rx(maxca,maxca)
3089       real*8 ry(maxca,maxca),rz(maxca,maxca)
3090       real*8 bx(maxca),by(maxca),bz(maxca)
3091       real*8 dis(maxca,maxca)
3092       real*8 shefx(maxca,12),shefy(maxca,12)
3093       real*8 shefz(maxca,12)
3094       real*8 dp45,dm45,w_beta
3095       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3096      $     c00,s00,ulnex,dnex
3097       integer    inb,nmax,iselect
3098 cc**********************************************************************
3099       common /phys1/     inb,nmax,iselect
3100       common /kyori2/    dis
3101       common /difvec/   rx,ry,rz
3102       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3103      $     c00,s00,ulnex,dnex
3104       common /sheconst/ dp45,dm45,w_beta
3105       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3106       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3107       common /shef/   shefx,shefy,shefz
3108 ci      integer istrand(maxca,maxca)
3109 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3110 ci      common  /shetest/ istrand,istrand_p,istrand_m
3111 cc**********************************************************************
3112 C     local variables
3113       integer  i,imm,im,jp,jpp,j,ip
3114       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3115       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3116       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3117       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3118       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
3119       real*8  yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
3120       real*8 e_gcont,fprim_gcont
3121 C********************************************************************************
3122       do i=2,inb-6
3123          ip=i+1
3124          im=i-1
3125          do j=i+3,inb-3
3126
3127             if (dis(im,j).lt.dfa_cutoff) then
3128             call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3129      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3130
3131             jp=j+1
3132             jpp=j+2
3133
3134 ci        if(istrand(im,j).eq.1
3135 ci     &    .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
3136
3137
3138             yy1=-(dis(i,jp)-ulhb)/dlhb
3139             y1x=rx(jp,i)/dis(i,jp)
3140             y1y=ry(jp,i)/dis(i,jp)
3141             y1z=rz(jp,i)/dis(i,jp)
3142             y11x=yy1*y1x
3143             y11y=yy1*y1y
3144             y11z=yy1*y1z
3145
3146             yy33=1.0D0/(dis(i,jp)*dis(i,ip))
3147             yyy3a=pin1(im,j)/(dis(i,jp)**2)
3148             yyy3b=pin1(im,j)/(dis(i,ip)**2)
3149             yy3=-pin1(im,j)/dshe
3150             y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
3151      $           +yyy3b*rx(i,ip))*yy3
3152             y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
3153      $           +yyy3b*ry(i,ip))*yy3
3154             y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
3155      $           +yyy3b*rz(i,ip))*yy3
3156
3157             yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3158             yyy4=pin2(im,j)/(dis(i,jp)**2)
3159             yy4=-pin2(im,j)/dshe
3160             y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3161             y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3162             y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3163
3164             yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3165             yyy5=pin3(im,j)/(dis(i,ip)**2)
3166             yy5=-pin3(im,j)/dshe
3167             y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3168             y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3169             y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3170
3171             sx=y11x+y3x+y4x+y5x
3172             sy=y11y+y3y+y4y+y5y
3173             sz=y11z+y3z+y4z+y5z
3174
3175             sx1=y11x+y3x+y4x
3176             sy1=y11y+y3y+y4y
3177             sz1=y11z+y3z+y4z
3178             sx2=y5x
3179             sy2=y5y
3180             sz2=y5z
3181
3182             shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3183      $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3184             shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3185      $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3186             shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3187      $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3188 !            shefx(i,6)=shefx(i,6)
3189 !     $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3190 !            shefy(i,6)=shefy(i,6)
3191 !     $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3192 !            shefz(i,6)=shefz(i,6)
3193 !     $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3194
3195             yy6=-(dis(jpp,i)-uldhb)/dldhb
3196             y6x=rx(jpp,i)/dis(jpp,i)
3197             y6y=ry(jpp,i)/dis(jpp,i)
3198             y6z=rz(jpp,i)/dis(jpp,i)
3199             y66x=yy6*y6x
3200             y66y=yy6*y6y
3201             y66z=yy6*y6z
3202
3203             yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3204             yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3205             yyy8b=pina1(im,j)/(dis(i,ip)**2)
3206             yy8=-pina1(im,j)/dshe
3207             y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3208      $           +yyy8b*rx(i,ip))*yy8
3209             y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3210      $           +yyy8b*ry(i,ip))*yy8
3211             y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3212      $           +yyy8b*rz(i,ip))*yy8
3213
3214             yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3215             yyy9=pina2(im,j)/(dis(i,jpp)**2)
3216             yy9=-pina2(im,j)/dshe
3217             y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3218             y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3219             y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3220
3221             yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3222             yyy10=pina3(im,j)/(dis(i,ip)**2)
3223             yy10=-pina3(im,j)/dshe
3224             y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3225             y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3226             y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3227
3228             sx=y66x+y8x+y9x+y10x
3229             sy=y66y+y8y+y9y+y10y
3230             sz=y66z+y8z+y9z+y10z
3231
3232             sx1=y66x+y8x+y9x
3233             sy1=y66y+y8y+y9y
3234             sz1=y66z+y8z+y9z
3235             sx2=y10x
3236             sy2=y10y
3237             sz2=y10z
3238
3239             shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3240      $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3241            shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3242      $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3243             shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3244      $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3245
3246 !            shefx(i,6)=shefx(i,6)
3247 !     $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3248 !           shefy(i,6)=shefy(i,6)
3249 !     $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3250 !            shefz(i,6)=shefz(i,6)
3251 !     $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3252
3253             endif
3254 ci         endif
3255
3256          enddo
3257       enddo
3258
3259       return
3260       end
3261 c-----------------------------------------------------------------------
3262       subroutine sheetforce11
3263       implicit none
3264       integer maxca
3265       parameter(maxca=800)
3266       real*8 dfa_cutoff,dfa_cutoff_delta
3267       parameter(dfa_cutoff=15.5d0)
3268       parameter(dfa_cutoff_delta=0.5d0)
3269 cc**********************************************************************
3270       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3271       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3272       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3273       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3274       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3275       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3276       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3277       real*8 rx(maxca,maxca)
3278       real*8 ry(maxca,maxca),rz(maxca,maxca)
3279       real*8 bx(maxca),by(maxca),bz(maxca)
3280       real*8 dis(maxca,maxca)
3281       real*8 shefx(maxca,12),shefy(maxca,12)
3282       real*8 shefz(maxca,12)
3283       real*8 dp45,dm45,w_beta
3284       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3285      $     c00,s00,ulnex,dnex
3286       integer    inb,nmax,iselect
3287 cc**********************************************************************
3288       common /phys1/     inb,nmax,iselect
3289       common /kyori2/    dis
3290       common /difvec/   rx,ry,rz
3291       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3292      $     c00,s00,ulnex,dnex
3293       common /sheconst/ dp45,dm45,w_beta
3294       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3295       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3296       common /shef/   shefx,shefy,shefz
3297 ci      integer istrand(maxca,maxca)
3298 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3299 ci      common  /shetest/ istrand,istrand_p,istrand_m
3300 C********************************************************************************
3301 C     local variables
3302       integer  j,jm,jmm,ip,i,ipp
3303       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3304       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3305       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3306       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3307       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3308       real*8  yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3309       real*8 e_gcont,fprim_gcont
3310 C********************************************************************************
3311
3312       do j=7,inb-1
3313          jm=j-1
3314          jmm=j-2
3315          do i=1,j-6
3316
3317             if (dis(i,jmm).lt.dfa_cutoff) then
3318             call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3319      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3320
3321             ip=i+1
3322             ipp=i+2
3323
3324 ci            if(istrand(i,jmm).eq.1
3325 ci     &   .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3326
3327
3328             yy1=-(dis(ipp,j)-ulhb)/dlhb
3329             y1x=rx(ipp,j)/dis(ipp,j)
3330             y1y=ry(ipp,j)/dis(ipp,j)
3331             y1z=rz(ipp,j)/dis(ipp,j)
3332             y11x=yy1*y1x
3333             y11y=yy1*y1y
3334             y11z=yy1*y1z
3335
3336             yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3337             yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3338             yy3=-pin2(i,jmm)/dshe
3339             y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3340             y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3341             y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3342
3343             yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3344             yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3345             yy4=-pin3(i,jmm)/dshe
3346             y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3347             y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3348             y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3349
3350             yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3351             yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3352             yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3353             yy5=-pin4(i,jmm)/dshe
3354             y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3355      $           -yyy5b*rx(jm,j))*yy5
3356             y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3357      $           -yyy5b*ry(jm,j))*yy5
3358             y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3359      $           -yyy5b*rz(jm,j))*yy5
3360
3361             sx=y11x+y3x+y4x+y5x
3362             sy=y11y+y3y+y4y+y5y
3363             sz=y11z+y3z+y4z+y5z
3364
3365             sx1=y3x
3366             sy1=y3y
3367             sz1=y3z
3368             sx2=y11x+y4x+y5x
3369             sy2=y11y+y4y+y5y
3370             sz2=y11z+y4z+y5z
3371
3372             shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3373      $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3374             shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3375      $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3376             shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3377      $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3378
3379 !            shefx(j,11)=shefx(j,11)
3380 !     $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3381 !            shefy(j,11)=shefy(j,11)
3382 !     $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3383 !            shefz(j,11)=shefz(j,11)
3384 !     $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3385
3386             yy6=-(dis(ip,j)-uldhb)/dldhb
3387             y6x=rx(ip,j)/dis(ip,j)
3388             y6y=ry(ip,j)/dis(ip,j)
3389             y6z=rz(ip,j)/dis(ip,j)
3390             y66x=yy6*y6x
3391             y66y=yy6*y6y
3392             y66z=yy6*y6z
3393
3394             yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3395             yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3396             yy8=-pina1(i,jmm)/dshe
3397             y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3398             y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3399             y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3400
3401             yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3402             yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3403             yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3404             yy9=-pina2(i,jmm)/dshe
3405             y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3406      $           -yyy9b*rx(jm,j))*yy9
3407             y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3408      $           -yyy9b*ry(jm,j))*yy9
3409             y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3410      $           -yyy9b*rz(jm,j))*yy9
3411
3412             yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3413             yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3414             yy10=-pina4(i,jmm)/dshe
3415             y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3416             y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3417             y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3418
3419             sx=y66x+y8x+y9x+y10x
3420             sy=y66y+y8y+y9y+y10y
3421             sz=y66z+y8z+y9z+y10z
3422
3423             sx1=y66x+y8x+y9x
3424             sy1=y66y+y8y+y9y
3425             sz1=y66z+y8z+y9z
3426             sx2=y10x
3427             sy2=y10y
3428             sz2=y10z
3429
3430             shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3431      $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3432            shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3433      $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3434             shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3435      $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3436
3437 !            shefx(j,11)=shefx(j,11)
3438 !     $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3439 !            shefy(j,11)=shefy(j,11)
3440 !     $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3441 !            shefz(j,11)=shefz(j,11)
3442 !     $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3443
3444             endif
3445 ci         endif
3446
3447          enddo
3448       enddo
3449
3450       return
3451       end
3452 c-----------------------------------------------------------------------
3453       subroutine sheetforce12
3454       implicit none
3455       integer maxca
3456       parameter(maxca=800)
3457       real*8 dfa_cutoff,dfa_cutoff_delta
3458       parameter(dfa_cutoff=15.5d0)
3459       parameter(dfa_cutoff_delta=0.5d0)
3460 cc**********************************************************************
3461       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3462       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3463       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3464       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3465       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3466       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3467       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3468       real*8 rx(maxca,maxca)
3469       real*8 ry(maxca,maxca),rz(maxca,maxca)
3470       real*8 bx(maxca),by(maxca),bz(maxca)
3471       real*8 dis(maxca,maxca)
3472       real*8 shefx(maxca,12),shefy(maxca,12)
3473       real*8 shefz(maxca,12)
3474       real*8 dp45,dm45,w_beta
3475       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3476      $     c00,s00,ulnex,dnex
3477       integer    inb,nmax,iselect
3478 cc**********************************************************************
3479       common /phys1/     inb,nmax,iselect
3480       common /kyori2/    dis
3481       common /difvec/   rx,ry,rz
3482       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3483      $     c00,s00,ulnex,dnex
3484       common /sheconst/ dp45,dm45,w_beta
3485       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3486       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3487       common /shef/   shefx,shefy,shefz
3488 ci      integer istrand(maxca,maxca)
3489 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3490 ci      common  /shetest/ istrand,istrand_p,istrand_m
3491 cc**********************************************************************
3492 C     local variables
3493       integer j,jm,jmm,ip,i,ipp,jp
3494       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3495       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3496       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3497       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3498       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3499       real*8 e_gcont,fprim_gcont
3500 !c*************************************************************************c
3501       do j=6,inb-2
3502          jp=j+1
3503          jm=j-1
3504          do i=1,j-5
3505
3506             if (dis(i,jm).lt.dfa_cutoff) then
3507             call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
3508      &                     dfa_cutoff_delta,e_gcont,fprim_gcont)
3509
3510             ip=i+1
3511             ipp=i+2
3512
3513 ci            if(istrand(i,jm).eq.1
3514 ci     &   .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3515
3516
3517             yy1=-(dis(ip,j)-ulhb)/dlhb
3518             y1x=rx(ip,j)/dis(ip,j)
3519             y1y=ry(ip,j)/dis(ip,j)
3520             y1z=rz(ip,j)/dis(ip,j)
3521             y11x=y1x*yy1
3522             y11y=y1y*yy1
3523             y11z=y1z*yy1
3524
3525             yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3526             yyy3=pin1(i,jm)/(dis(ip,j)**2)
3527             yy3=-pin1(i,jm)/dshe
3528             y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3529             y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3530             y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3531             yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3532
3533             yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3534             yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3535             yy4=-pin2(i,jm)/dshe
3536             y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3537      $           +yyy4b*rx(j,jp))*yy4
3538             y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3539      $           +yyy4b*ry(j,jp))*yy4
3540             y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3541      $           +yyy4b*rz(j,jp))*yy4
3542
3543             yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3544             yyy5=pin4(i,jm)/(dis(j,jp)**2)
3545             yy5=-pin4(i,jm)/dshe
3546             y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3547             y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3548             y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3549
3550             sx=y11x+y3x+y4x+y5x
3551             sy=y11y+y3y+y4y+y5y
3552             sz=y11z+y3z+y4z+y5z
3553
3554             sx1=y11x+y3x+y4x
3555             sy1=y11y+y3y+y4y
3556             sz1=y11z+y3z+y4z
3557             sx2=y5x
3558             sy2=y5y
3559             sz2=y5z
3560
3561             shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3562      $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3563             shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3564      $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3565             shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3566      $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3567
3568 !            shefx(j,12)=shefx(j,12)
3569 !     $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3570 !            shefy(j,12)=shefy(j,12)
3571 !     $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3572 !            shefz(j,12)=shefz(j,12)
3573 !     $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3574
3575             yy6=-(dis(ipp,j)-uldhb)/dldhb
3576             y6x=rx(ipp,j)/dis(ipp,j)
3577             y6y=ry(ipp,j)/dis(ipp,j)
3578             y6z=rz(ipp,j)/dis(ipp,j)
3579             y66x=yy6*y6x
3580             y66y=yy6*y6y
3581             y66z=yy6*y6z
3582
3583             yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3584             yyy8=pina2(i,jm)/(dis(j,jp)**2)
3585             yy8=-pina2(i,jm)/dshe
3586             y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3587             y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3588             y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3589
3590             yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3591             yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3592             yy9=-pina3(i,jm)/dshe
3593             y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3594             y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3595             y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3596
3597             yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3598             yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3599             yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3600             yy10=-pina4(i,jm)/dshe
3601             y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3602      $           +yyy10b*rx(j,jp))*yy10
3603             y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3604      $           +yyy10b*ry(j,jp))*yy10
3605             y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3606      $           +yyy10b*rz(j,jp))*yy10
3607
3608             sx=y66x+y8x+y9x+y10x
3609             sy=y66y+y8y+y9y+y10y
3610             sz=y66z+y8z+y9z+y10z
3611
3612             sx1=y8x
3613             sy1=y8y
3614             sz1=y8z
3615             sx2=y66x+y9x+y10x
3616             sy2=y66y+y9y+y10y
3617             sz2=y66z+y9z+y10z
3618
3619             shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3620      $           -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3621            shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3622      $           -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3623             shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3624      $           -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3625
3626             endif
3627
3628 ci         endif
3629
3630          ENDDO
3631       ENDDO
3632
3633       RETURN
3634       END
3635 C===============================================================================