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