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