source/wham/src/dfa.F

   1       subroutine init_dfa_vars
   2
   3       include 'DIMENSIONS'
   4       include 'COMMON.INTERACT'
   5       include 'COMMON.DFA'
   6
   7       integer ii
   8
   9 C     Number of restraints
  10       idisnum = 0
  11       iphinum = 0
  12       ithenum = 0
  13       ineinum = 0
  14
  15       idislis = 0
  16       iphilis = 0
  17       ithelis = 0
  18       ineilis = 0
  19       jneilis = 0
  20       jneinum = 0
  21       kshell  = 0
  22       fnei    = 0
  23 C     For beta
  24       nca     = 0
  25       icaidx  = 0
  26
  27 C     real variables
  28 CC    WEIGHTS for each min
  29       sccdist = 0.0d0
  30       fdist   = 0.0d0
  31       sccphi  = 0.0d0
  32       sccthe  = 0.0d0
  33       sccnei  = 0.0d0
  34       fphi1   = 0.0d0
  35       fphi2   = 0.0d0
  36       fthe1   = 0.0d0
  37       fthe2   = 0.0d0
  38 C     energies
  39       edfatot = 0.0d0
  40       edfadis = 0.0d0
  41       edfaphi = 0.0d0
  42       edfathe = 0.0d0
  43       edfanei = 0.0d0
  44       edfabet = 0.0d0
  45 C     weights for each E term
  46 C     these should be identical with
  47       dis_inc = 0.0d0
  48       phi_inc = 0.0d0
  49       the_inc = 0.0d0
  50       nei_inc = 0.0d0
  51       beta_inc = 0.0d0
  52       wshet   = 0.0d0
  53 C     precalculate exp table!
  54 c      dfaexp  = 0.0d0
  55 c      do ii = 1, 15001
  56 c         dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
  57 c      end do
  58
  59       ishiftca=nnt-1
  60       ilastca=nct
  61
  62       print *,'ishiftca=',ishiftca,'ilastca=',ilastca
  63
  64       return
  65       end
  66
  67
  68       subroutine read_dfa_info
  69 C
  70 C     read fragment informations
  71 C
  72       implicit real*8 (a-h,o-z)
  73       include 'DIMENSIONS'
  74       include '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       vbeta=vbeta*beta_inc
1089       enebet=vbeta
1090       edfabeta=enebet
1091
1092       do j=1,nca
1093          gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1094          gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1095          gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1096       enddo
1097
1098 #ifdef DEBUG1
1099       do j=1,nca
1100         write(*,'(5x,i5,10x,3f10.5)') j,bx(j),by(j),bz(j)
1101       enddo
1102
1103
1104       gdfab=0
1105       dinc=0.001
1106       do j=1,nca
1107         cth_orig(j)=cth(j)
1108         sth_orig(j)=sth(j)
1109       enddo
1110
1111       do j=1,nca
1112
1113        bx(j)=bx(j)+dinc
1114        call angvectors(nca)
1115        bx(j)=bx(j)-2*dinc
1116        call angvectors(nca)
1117        atxnum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1118        astxnum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1119        if (j.gt.1) then
1120        atmxnum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1121        astmxnum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1122        endif
1123        if (j.gt.2) then
1124        atmmxnum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1125        astmmxnum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1126        endif
1127        if (j.gt.3) then
1128        atm3xnum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1129        astm3xnum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1130        endif
1131        bx(j)=bx(j)+dinc
1132        by(j)=by(j)+dinc
1133        call angvectors(nca)
1134        by(j)=by(j)-2*dinc
1135        call angvectors(nca)
1136        by(j)=by(j)+dinc
1137        atynum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1138        astynum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1139        if (j.gt.1) then
1140        atmynum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1141        astmynum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1142        endif
1143        if (j.gt.2) then
1144        atmmynum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1145        astmmynum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1146        endif
1147        if (j.gt.3) then
1148        atm3ynum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1149        astm3ynum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1150        endif
1151
1152        bz(j)=bz(j)+dinc
1153        call angvectors(nca)
1154        bz(j)=bz(j)-2*dinc
1155        call angvectors(nca)
1156        bz(j)=bz(j)+dinc
1157
1158        atznum(j)=0.5*(cth(j)-cth_orig(j))/dinc
1159        astznum(j)=0.5*(sth(j)-sth_orig(j))/dinc
1160        if (j.gt.1) then
1161        atmznum(j)=0.5*(cth(j-1)-cth_orig(j-1))/dinc
1162        astmznum(j)=0.5*(sth(j-1)-sth_orig(j-1))/dinc
1163        endif
1164        if (j.gt.2) then
1165        atmmznum(j)=0.5*(cth(j-2)-cth_orig(j-2))/dinc
1166        astmmznum(j)=0.5*(sth(j-2)-sth_orig(j-2))/dinc
1167        endif
1168        if (j.gt.3) then
1169        atm3znum(j)=0.5*(cth(j-3)-cth_orig(j-3))/dinc
1170        astm3znum(j)=0.5*(sth(j-3)-sth_orig(j-3))/dinc
1171        endif
1172
1173       enddo
1174
1175       do i=1,nca
1176         write (*,'(2i5,a2,6f10.5)')
1177      &  i,1,"x",atxnum(i),atx(i),atxnum(i)/atx(i),
1178      &          astxnum(i),astx(i),astxnum(i)/astx(i),
1179      &  i,1,"y",atynum(i),aty(i),atynum(i)/aty(i),
1180      &          astynum(i),asty(i),astynum(i)/asty(i),
1181      &  i,1,"z",atznum(i),atz(i),atznum(i)/atz(i),
1182      &          astznum(i),astz(i),astznum(i)/astz(i),
1183      &  i,2,"x",atmxnum(i),atmx(i),atmxnum(i)/atmx(i),
1184      &          astmxnum(i),astmx(i),astmxnum(i)/astmx(i),
1185      &  i,2,"y",atmynum(i),atmy(i),atmynum(i)/atmy(i),
1186      &          astmynum(i),astmy(i),astmynum(i)/astmy(i),
1187      &  i,2,"z",atmznum(i),atmz(i),atmznum(i)/atmz(i),
1188      &          astmznum(i),astmz(i),astmznum(i)/astmz(i),
1189      &  i,3,"x",atmmxnum(i),atmmx(i),atmmxnum(i)/atmmx(i),
1190      &          astmmxnum(i),astmmx(i),astmmxnum(i)/astmmx(i),
1191      &  i,3,"y",atmmynum(i),atmmy(i),atmmynum(i)/atmmy(i),
1192      &          astmmynum(i),astmmy(i),astmmynum(i)/astmmy(i),
1193      &  i,3,"z",atmmznum(i),atmmz(i),atmmznum(i)/atmmz(i),
1194      &          astmmznum(i),astmmz(i),astmmznum(i)/astmmz(i),
1195      &  i,4,"x",atm3xnum(i),atm3x(i),atm3xnum(i)/atm3x(i),
1196      &          astm3xnum(i),astm3x(i),astm3xnum(i)/astm3x(i),
1197      &  i,4,"y",atm3ynum(i),atm3y(i),atm3ynum(i)/atm3y(i),
1198      &          astm3ynum(i),astm3y(i),astm3ynum(i)/astm3y(i),
1199      &  i,4,"z",atm3znum(i),atm3z(i),atm3znum(i)/atm3z(i),
1200      &          astm3znum(i),astm3z(i),astm3znum(i)/astm3z(i),
1201      &  i,0," ",cth_orig(i),sth_orig(i)
1202       enddo
1203
1204
1205       gdfab=0
1206       dinc=0.001
1207
1208       do j=1,nca
1209
1210        bx(j)=bx(j)+dinc
1211        call angvectors(nca)
1212        call sheetforce(nca,wshet)
1213        vbeta1=vbeta*beta_inc
1214        bx(j)=bx(j)-2*dinc
1215        call angvectors(nca)
1216        call sheetforce(nca,wshet)
1217        vbeta2=vbeta*beta_inc
1218        gdfab(1,j)=(vbeta2-vbeta1)/dinc/2
1219        bx(j)=bx(j)+dinc
1220
1221        by(j)=by(j)+dinc
1222        call angvectors(nca)
1223        call sheetforce(nca,wshet)
1224        vbeta1=vbeta*beta_inc
1225        by(j)=by(j)-2*dinc
1226        call angvectors(nca)
1227        call sheetforce(nca,wshet)
1228        vbeta2=vbeta*beta_inc
1229        gdfab(2,j)=(vbeta2-vbeta1)/dinc/2
1230        by(j)=by(j)+dinc
1231
1232        bz(j)=bz(j)+dinc
1233        call angvectors(nca)
1234        call sheetforce(nca,wshet)
1235        vbeta1=vbeta*beta_inc
1236        bz(j)=bz(j)-2*dinc
1237        call angvectors(nca)
1238        call sheetforce(nca,wshet)
1239        vbeta2=vbeta*beta_inc
1240        gdfab(3,j)=(vbeta2-vbeta1)/dinc/2
1241        bz(j)=bz(j)+dinc
1242
1243
1244       enddo
1245
1246
1247       call angvectors(nca)
1248       call sheetforce(nca,wshet)
1249       do j=1,nca
1250          shetfx(j)=shetfx(j)*beta_inc
1251          shetfy(j)=shetfy(j)*beta_inc
1252          shetfz(j)=shetfz(j)*beta_inc
1253       enddo
1254
1255
1256       write(*,*) 'xyz analytical and numerical gradient'
1257       do j=1,nca
1258         write(*,'(5x,i5,10x,6f10.5)') j,-shetfx(j),-shetfy(j),-shetfz(j)
1259      &                   ,(-gdfab(i,j),i=1,3)
1260       enddo
1261
1262       do j=1,nca
1263         write(*,'(5x,i5,10x,3f10.2)') j,shetfx(j)/gdfab(1,j),
1264      &                                  shetfy(j)/gdfab(2,j),
1265      &                                  shetfz(j)/gdfab(3,j)
1266       enddo
1267
1268       stop
1269 #endif
1270
1271       return
1272       end
1273 C-------------------------------------------------------------------------------
1274       subroutine angvectors(nca)
1275 c      implicit real*4(a-h,o-z)
1276       implicit none
1277       integer nca
1278       integer maxca
1279       parameter(maxca=800)
1280       real*8   pai,zero
1281       parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1282
1283       real*8   bx(maxca),by(maxca),bz(maxca)
1284       real*8   dis(maxca,maxca)
1285       real*8   apx(maxca),apy(maxca),apz(maxca)
1286       real*8   apmx(maxca),apmy(maxca),apmz(maxca)
1287       real*8   apmmx(maxca),apmmy(maxca),apmmz(maxca)
1288       real*8   apm3x(maxca),apm3y(maxca),apm3z(maxca)
1289       real*8   atx(maxca),aty(maxca),atz(maxca)
1290       real*8   atmx(maxca),atmy(maxca),atmz(maxca)
1291       real*8   atmmx(maxca),atmmy(maxca),atmmz(maxca)
1292       real*8   atm3x(maxca),atm3y(maxca),atm3z(maxca)
1293       real*8   astx(maxca),asty(maxca),astz(maxca)
1294       real*8   astmx(maxca),astmy(maxca),astmz(maxca)
1295       real*8   astmmx(maxca),astmmy(maxca),astmmz(maxca)
1296       real*8   astm3x(maxca),astm3y(maxca),astm3z(maxca)
1297       real*8   sth(maxca)
1298       real*8   cph(maxca),cth(maxca)
1299       real*8   ulcos(maxca)
1300       real*8   p,c
1301       integer  i, ip, ipp, ip3, j
1302       real*8   rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1303       real*8   rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1304       real*8   gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1305       real*8   cix, ciy, ciz, cipx, cipy, cipz
1306       real*8   gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1307       real*8   d10, d11, d12, d13, d20, d21, d22, d23, d24
1308       real*8   d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1309       real*8   d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1310       real*8   dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1311       real*8   dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1312       real*8   g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1313       real*8   gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1314       real*8   gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1315       real*8   gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1316       real*8   gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1317       real*8   grpp,gx,gy,gz
1318       real*8   rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1319       real*8   sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1320       integer inb,nmax,iselect
1321
1322       common /sheca/   bx,by,bz
1323       common /difvec/  rx, ry, rz
1324       common /ulang/    ulcos
1325       common /phys1/   inb,nmax,iselect
1326       common /phys4/   p,c
1327       common /kyori2/  dis
1328       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1329      &     apmmz,apm3x,apm3y,apm3z
1330       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1331      &     atmmz,atm3x,atm3y,atm3z
1332       common /coscos/   cph,cth
1333       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1334      &     astmmz,astm3x,astm3y,astm3z
1335       common /sinsin/   sth
1336 C-------------------------------------------------------------------------------
1337 c      write(*,*) 'inside angvectors'
1338 C     initialize
1339       p=0.1d0
1340       c=1.0d0
1341       inb=nca
1342       cph=zero; cth=zero; sth=zero
1343       apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1344       apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1345       atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1346       atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1347       astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1348       astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1349       astm3z=zero
1350 C     end of initialize
1351 C     r[x,y,z] calc and distance calculation
1352       rx=zero;ry=zero;rz=zero
1353
1354       do i=1,inb
1355          do j=1,inb
1356             rx(i,j)=bx(j)-bx(i)
1357             ry(i,j)=by(j)-by(i)
1358             rz(i,j)=bz(j)-bz(i)
1359             dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1360 c            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1361 c            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1362 c            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1363 c            write(*,*) 'dis(i,j):',i,j,dis(i,j)
1364          enddo
1365       enddo
1366 c     end of r[x,y,z] calc
1367 C     cos calc
1368       do i=1,inb-2
1369          ip=i+1
1370          ipp=i+2
1371
1372          if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1373             ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1374      $           +rz(i,ip)*rz(ip,ipp)
1375             ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1376          endif
1377       enddo
1378 c     end of virtual bond angle
1379 c      write(*,*) 'inside angvectors1'
1380 crc       do i=1,inb-3
1381       do i=1,inb
1382          ip=i+1
1383          ipp=i+2
1384          ip3=i+3
1385          rix=bx(ip)-bx(i)
1386          riy=by(ip)-by(i)
1387          riz=bz(ip)-bz(i)
1388          ripx=bx(ipp)-bx(ip)
1389          ripy=by(ipp)-by(ip)
1390          ripz=bz(ipp)-bz(ip)
1391          rippx=bx(ip3)-bx(ipp)
1392          rippy=by(ip3)-by(ipp)
1393          rippz=bz(ip3)-bz(ipp)
1394
1395          gx=riy*ripz-riz*ripy
1396          gy=riz*ripx-rix*ripz
1397          gz=rix*ripy-riy*ripx
1398          gpx=ripy*rippz-ripz*rippy
1399          gpy=ripz*rippx-ripx*rippz
1400          gpz=ripx*rippy-ripy*rippx
1401          gpcrp_x=gpy*ripz-gpz*ripy
1402          gpcrp_y=gpz*ripx-gpx*ripz
1403          gpcrp_z=gpx*ripy-gpy*ripx
1404          d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1405          gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1406      &        -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1407
1408          if(i.ge.2) then
1409             rimx=bx(i)-bx(i-1)
1410             rimy=by(i)-by(i-1)
1411             rimz=bz(i)-bz(i-1)
1412             gmx=rimy*riz-rimz*riy
1413             gmy=rimz*rix-rimx*riz
1414             gmz=rimx*riy-rimy*rix
1415             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1416             dgm3=dgm**3
1417             ggm=gmx*gx+gmy*gy+gmz*gz
1418             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1419             drim=dis(i-1,i)
1420             drim3=drim**3
1421             gcr_x=gy*riz-gz*riy
1422             gcr_y=gz*rix-gx*riz
1423             gcr_z=gx*riy-gy*rix
1424             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1425             d_gcr3=d_gcr**3
1426             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1427      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1428          endif
1429 c         write(*,*) 'inside angvectors2'
1430          if(i.ge.3) then
1431             rimmx=bx(i-1)-bx(i-2)
1432             rimmy=by(i-1)-by(i-2)
1433             rimmz=bz(i-1)-bz(i-2)
1434             drimm=dis(i-2,i-1)
1435             gmmx=rimmy*rimz-rimmz*rimy
1436             gmmy=rimmz*rimx-rimmx*rimz
1437             gmmz=rimmx*rimy-rimmy*rimx
1438             dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1439             dgmm3=dgmm**3
1440             gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1441             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1442             gmcrim_x=gmy*rimz-gmz*rimy
1443             gmcrim_y=gmz*rimx-gmx*rimz
1444             gmcrim_z=gmx*rimy-gmy*rimx
1445             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1446             d_gmcrim3=d_gmcrim**3
1447             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1448      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1449          endif
1450
1451          if(i.ge.4) then
1452             rim3x=bx(i-2)-bx(i-3)
1453             rim3y=by(i-2)-by(i-3)
1454             rim3z=bz(i-2)-bz(i-3)
1455             g3x=rim3y*rimmz-rim3z*rimmy
1456             g3y=rim3z*rimmx-rim3x*rimmz
1457             g3z=rim3x*rimmy-rim3y*rimmx
1458             dg30=sqrt(g3x**2+g3y**2+g3z**2)
1459             g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1460             g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1461 cc**********************************************************************
1462             gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1463             gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1464             gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1465             d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1466             d_gmmcrimm3=d_gmmcrimm**3
1467             gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1468      &           -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1469          endif
1470
1471          dri=dis(i,i+1)
1472          drip=dis(i+1,i+2)
1473          dripp=dis(i+2,i+3)
1474          dri3=dri**3
1475          dg=sqrt(gx**2+gy**2+gz**2)
1476          dgp=sqrt(gpx**2+gpy**2+gpz**2)
1477          dg3=dg**3
1478
1479          ggp=gx*gpx+gy*gpy+gz*gpz
1480          grpp=gx*rippx+gy*rippy+gz*rippz
1481
1482          if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1483      &        .and.d_gpcrp.gt.0.0D0) then
1484             cph(i)=grpp/dg/dripp
1485             cth(i)=ggp/dg/dgp
1486             sth(i)=gpcrp__g/d_gpcrp/dg
1487          else
1488 c
1489             cph(i)=1.0D0
1490             cth(i)=1.0D0
1491             sth(i)=0.0D0
1492          endif
1493
1494 c         write(*,*) 'inside angvectors3'
1495
1496          if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1497      &        .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1498             d10=1.0D0/(dg*dgp)
1499             d11=ggp/(dg3*dgp)
1500             d12=1.0D0/(dg*dripp)
1501             d13=grpp/(dg3*dripp)
1502             sd10=1.0D0/(d_gpcrp*dg)
1503             sd11=gpcrp__g/(d_gpcrp*dg3)
1504          else
1505             d10=0.0D0
1506             d11=0.0D0
1507             d12=0.0D0
1508             d13=0.0D0
1509             sd10=0.0D0
1510             sd11=0.0D0
1511          endif
1512
1513          atx(i)=(ripz*gpy-ripy*gpz)*d10
1514      &        -(gy*ripz-gz*ripy)*d11
1515          aty(i)=(ripx*gpz-ripz*gpx)*d10
1516      &        -(gz*ripx-gx*ripz)*d11
1517          atz(i)=(ripy*gpx-ripx*gpy)*d10
1518      &        -(gx*ripy-gy*ripx)*d11
1519          astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1520      &        +ripy*gpy*ripx-gpx*ripz**2)
1521      &        -sd11*(gy*ripz-gz*ripy)
1522          asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1523      &        -gpy*ripx**2+gpz*ripy*ripz)
1524      &        -sd11*(-gx*ripz+gz*ripx)
1525          astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1526      &        -gpz*ripy**2+ripz*gpx*ripx)
1527      &        -sd11*(gx*ripy-gy*ripx)
1528          apx(i)=(ripz*rippy-ripy*rippz)*d12
1529      &        -(gy*ripz-gz*ripy)*d13
1530          apy(i)=(ripx*rippz-ripz*rippx)*d12
1531      &        -(gz*ripx-gx*ripz)*d13
1532          apz(i)=(ripy*rippx-ripx*rippy)*d12
1533      &        -(gx*ripy-gy*ripx)*d13
1534
1535          if(i.ge.2) then
1536             cix=bx(ip)-bx(i-1)
1537             ciy=by(ip)-by(i-1)
1538             ciz=bz(ip)-bz(i-1)
1539             cipx=bx(ipp)-bx(i)
1540             cipy=by(ipp)-by(i)
1541             cipz=bz(ipp)-bz(i)
1542             ripx=bx(ipp)-bx(ip)
1543             ripy=by(ipp)-by(ip)
1544             ripz=bz(ipp)-bz(ip)
1545             if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1546      &           .and.d_gcr3.gt.0.0D0) then
1547                d20=1.0D0/(dg*dgm)
1548                d21=ggm/(dgm3*dg)
1549                d22=ggm/(dgm*dg3)
1550                d23=1.0D0/(dgm*drip)
1551                d24=gmrp/(dgm3*drip)
1552                sd20=1.0D0/(d_gcr*dgm)
1553                sd21=gcr__gm/(d_gcr3*dgm)
1554                sd22=gcr__gm/(d_gcr*dgm3)
1555             else
1556                d20=0.0D0
1557                d21=0.0D0
1558                d22=0.0D0
1559                d23=0.0D0
1560                d24=0.0D0
1561                sd20=0.0D0
1562                sd21=0.0D0
1563                sd22=0.0D0
1564             endif
1565             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1566      &           -(ciy*gmz-ciz*gmy)*d21
1567      &           +(ripy*gz-ripz*gy)*d22
1568             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1569      &           -(ciz*gmx-cix*gmz)*d21
1570      &           +(ripz*gx-ripx*gz)*d22
1571             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1572      &           -(cix*gmy-ciy*gmx)*d21
1573      &           +(ripx*gy-ripy*gx)*d22
1574 cc**********************************************************************
1575             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1576      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1577      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1578      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1579      &           +gcr_z*(-ripz*rix+gy))
1580      &           -sd22*(-gmy*ciz+gmz*ciy)
1581
1582             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1583      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1584      &           +riz*ripz*gmy)
1585      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1586      &           -gcr_z*(ripz*riy+gx))
1587      &           -sd22*(gmx*ciz-gmz*cix)
1588
1589             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1590      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1591      &           -riz*gx*cix)
1592      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1593      &           +gcr_z*(ripy*riy+ripx*rix))
1594      &           -sd22*(-gmx*ciy+gmy*cix)
1595 cc**********************************************************************
1596             apmx(i)=(ciy*ripz-ripy*ciz)*d23
1597      &           -(ciy*gmz-ciz*gmy)*d24
1598             apmy(i)=(ciz*ripx-ripz*cix)*d23
1599      &           -(ciz*gmx-cix*gmz)*d24
1600             apmz(i)=(cix*ripy-ripx*ciy)*d23
1601      &           -(cix*gmy-ciy*gmx)*d24
1602          endif
1603
1604          if(i.ge.3) then
1605             if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1606      &           .and.d_gmcrim3.gt.0.0D0) then
1607                d30=1.0D0/(dgm*dgmm)
1608                d31=gmmgm/(dgm3*dgmm)
1609                d32=gmmgm/(dgm*dgmm3)
1610                d33=1.0D0/(dgmm*dri)
1611                d34=gmmr/(dgmm3*dri)
1612                d35=gmmr/(dgmm*dri3)
1613                sd30=1.0D0/(d_gmcrim*dgmm)
1614                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1615                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1616             else
1617                d30=0.0D0
1618                d31=0.0D0
1619                d32=0.0D0
1620                d33=0.0D0
1621                d34=0.0D0
1622                d35=0.0D0
1623                sd30=0.0D0
1624                sd31=0.0D0
1625                sd32=0.0D0
1626             endif
1627
1628 c            write(*,*) 'inside angvectors4'
1629
1630 cc**********************************************************************
1631             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1632      &           -(ciy*gmz-ciz*gmy)*d31
1633      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1634             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1635      &           -(ciz*gmx-cix*gmz)*d31
1636      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1637             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1638      &           -(cix*gmy-ciy*gmx)*d31
1639      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1640 cc**********************************************************************
1641             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1642      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1643      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1644      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1645      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1646      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1647      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1648
1649             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1650      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1651      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1652      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1653      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1654      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1655      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1656
1657             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1658      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1659      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1660      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1661      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1662      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1663      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1664 c**********************************************************************
1665             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1666      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1667      &           +rix*d35
1668             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1669      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1670      &           +riy*d35
1671             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1672      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1673      &           +riz*d35
1674          endif
1675
1676          if(i.ge.4) then
1677             if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1678      &           .and.drim3.gt.0.0D0
1679      &           .and.d_gmmcrimm3.gt.0.0D0) then
1680                d40=1.0D0/(dg30*dgmm)
1681                d41=g3gmm/(dg30*dgmm3)
1682                d42=1.0D0/(dg30*drim)
1683                d43=g3rim_/(dg30*drim3)
1684                sd40=1.0D0/(dg30*d_gmmcrimm)
1685                sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1686             else
1687                d40=0.0D0
1688                d41=0.0D0
1689                d42=0.0D0
1690                d43=0.0D0
1691                sd40=0.0D0
1692                sd41=0.0D0
1693             endif
1694             atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1695      &           -(gmmy*rimmz-gmmz*rimmy)*d41
1696             atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1697      &           -(gmmz*rimmx-gmmx*rimmz)*d41
1698             atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1699      &           -(gmmx*rimmy-gmmy*rimmx)*d41
1700 cc**********************************************************************
1701             astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1702      &           -g3z*rimmz*rimmx+rimmy**2*g3x)
1703      &           -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1704      &           -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1705
1706             astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1707      &           -rimmx*rimmy*g3x+rimmz**2*g3y)
1708      &           -sd41*(-gmmcrimm_x*rimmx*rimmy
1709      &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmy)
1710
1711 c     &           +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1712
1713             astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1714      &           +g3z*rimmx**2-rimmz*rimmy*g3y)
1715      &           -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1716      &           +gmmcrimm_z*(rimmy**2+rimmx**2))
1717 c**********************************************************************
1718             apm3x(i)=g3x*d42-rimx*d43
1719             apm3y(i)=g3y*d42-rimy*d43
1720             apm3z(i)=g3z*d42-rimz*d43
1721          endif
1722       enddo
1723 c*******************************************************************************
1724
1725 c      write(*,*) 'inside angvectors5'
1726
1727 c       do i=inb-2,inb
1728        do i=1,0
1729          rimx=bx(i)-bx(i-1)
1730          rimy=by(i)-by(i-1)
1731          rimz=bz(i)-bz(i-1)
1732          rimmx=bx(i-1)-bx(i-2)
1733          rimmy=by(i-1)-by(i-2)
1734          rimmz=bz(i-1)-bz(i-2)
1735          rim3x=bx(i-2)-bx(i-3)
1736          rim3y=by(i-2)-by(i-3)
1737          rim3z=bz(i-2)-bz(i-3)
1738          gmmx=rimmy*rimz-rimmz*rimy
1739          gmmy=rimmz*rimx-rimmx*rimz
1740          gmmz=rimmx*rimy-rimmy*rimx
1741          g3x=rim3y*rimmz-rim3z*rimmy
1742          g3y=rim3z*rimmx-rim3x*rimmz
1743          g3z=rim3x*rimmy-rim3y*rimmx
1744
1745          dg30=sqrt(g3x**2+g3y**2+g3z**2)
1746          g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1747          dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1748          dgmm3=dgmm**3
1749          drim=dis(i-1,i)
1750          drimm=dis(i-2,i-1)
1751          drim3=drim**3
1752          g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1753 cc**********************************************************************
1754          gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1755          gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1756          gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1757          d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1758          d_gmmcrimm3=d_gmmcrimm**3
1759          gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1760      &        -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1761
1762          if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1763      &        .and.drim3.gt.0.0D0
1764      &        .and.d_gmmcrimm3.gt.0.0D0) then
1765             d40=1.0D0/(dg30*dgmm)
1766             d41=g3gmm/(dg30*dgmm3)
1767             d42=1.0D0/(dg30*drim)
1768             d43=g3rim_/(dg30*drim3)
1769             sd40=1.0D0/(dg30*d_gmmcrimm)
1770             sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1771          else
1772             d40=0.0D0
1773             d41=0.0D0
1774             d42=0.0D0
1775             d43=0.0D0
1776             sd40=0.0D0
1777             sd41=0.0D0
1778          endif
1779          atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1780      &        -(gmmy*rimmz-gmmz*rimmy)*d41
1781          atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1782      &        -(gmmz*rimmx-gmmx*rimmz)*d41
1783          atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1784      &        -(gmmx*rimmy-gmmy*rimmx)*d41
1785 cc**********************************************************************
1786          astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1787      &        -g3z*rimmz*rimmx+rimmy**2*g3x)
1788      &        -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1789      &        -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1790
1791          astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1792      &        -rimmx*rimmy*g3x+rimmz**2*g3y)
1793      &        -sd41*(-gmmcrimm_x*rimmx*rimmy
1794      &        +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1795
1796          astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1797      &        +g3z*rimmx**2-rimmz*rimmy*g3y)
1798      &        -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1799      &        +gmmcrimm_z*(rimmy**2+rimmx**2))
1800 cc**********************************************************************
1801          apm3x(i)=g3x*d42-rimx*d43
1802          apm3y(i)=g3y*d42-rimy*d43
1803          apm3z(i)=g3z*d42-rimz*d43
1804
1805          if(i.le.inb-1) then
1806             ip=i+1
1807             rix=bx(ip)-bx(i)
1808             riy=by(ip)-by(i)
1809             riz=bz(ip)-bz(i)
1810             cix=bx(ip)-bx(i-1)
1811             ciy=by(ip)-by(i-1)
1812             ciz=bz(ip)-bz(i-1)
1813             gmx=rimy*riz-rimz*riy
1814             gmy=rimz*rix-rimx*riz
1815             gmz=rimx*riy-rimy*rix
1816             dgm=sqrt(gmx**2+gmy**2+gmz**2)
1817             dgm3=dgm**3
1818             dri=dis(i,i+1)
1819             dri3=dri**3
1820             gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1821             gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1822             gmcrim_x=gmy*rimz-gmz*rimy
1823             gmcrim_y=gmz*rimx-gmx*rimz
1824             gmcrim_z=gmx*rimy-gmy*rimx
1825             d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1826             d_gmcrim3=d_gmcrim**3
1827             gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1828      &           -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1829
1830             if(dgm3.gt.0.0D0.and.
1831      &           dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1832      &           .and.d_gmcrim3.gt.0.0D0) then
1833                d30=1.0D0/(dgm*dgmm)
1834                d31=gmmgm/(dgm3*dgmm)
1835                d32=gmmgm/(dgm*dgmm3)
1836                d33=1.0D0/(dgmm*dri)
1837                d34=gmmr/(dgmm3*dri)
1838                d35=gmmr/(dgmm*dri3)
1839                sd30=1.0D0/(d_gmcrim*dgmm)
1840                sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1841                sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1842
1843             else
1844                d30=0.0D0
1845                d31=0.0D0
1846                d32=0.0D0
1847                d33=0.0D0
1848                d34=0.0D0
1849                d35=0.0D0
1850                sd30=0.0D0
1851                sd31=0.0D0
1852                sd32=0.0D0
1853             endif
1854 cc**********************************************************************
1855             atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1856      &           -(ciy*gmz-ciz*gmy)*d31
1857      &           -(gmmy*rimmz-gmmz*rimmy)*d32
1858             atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1859      &           -(ciz*gmx-cix*gmz)*d31
1860      &           -(gmmz*rimmx-gmmx*rimmz)*d32
1861             atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1862      &           -(cix*gmy-ciy*gmx)*d31
1863      &           -(gmmx*rimmy-gmmy*rimmx)*d32
1864 cc**********************************************************************
1865             astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1866      &           +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1867      &           +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1868      &           -ciy*rimy*gmmx-rimz*gmx*rimmz)
1869      &           -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1870      &           +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1871      &           -sd32*(gmmy*rimmz-rimmy*gmmz)
1872
1873             astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1874      &           +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1875      &           -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1876      &           +gmz*rimy*rimmz-rimz*ciz*gmmy)
1877      &           -sd31*(gmcrim_x*(cix*rimy-gmz)
1878      &           +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1879      &           -sd32*(-gmmx*rimmz+rimmx*gmmz)
1880
1881             astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1882      &           +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1883      &           -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1884      &           +rimz*ciy*gmmy+rimz*gmx*rimmx)
1885      &           -sd31*(gmcrim_x*(cix*rimz+gmy)
1886      &           +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1887      &           -sd32*(gmmx*rimmy-rimmx*gmmy)
1888 cc**********************************************************************
1889             apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1890      &           -(gmmy*rimmz-gmmz*rimmy)*d34
1891      &           +rix*d35
1892             apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1893      &           -(gmmz*rimmx-gmmx*rimmz)*d34
1894      &           +riy*d35
1895             apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1896      &           -(gmmx*rimmy-gmmy*rimmx)*d34
1897      &           +riz*d35
1898          endif
1899
1900 c         write(*,*) 'inside angvectors6'
1901
1902          if(i.eq.inb-2) then
1903             ipp=i+2
1904             ripx=bx(ipp)-bx(ip)
1905             ripy=by(ipp)-by(ip)
1906             ripz=bz(ipp)-bz(ip)
1907             cipx=bx(ipp)-bx(i)
1908             cipy=by(ipp)-by(i)
1909             cipz=bz(ipp)-bz(i)
1910             gx=riy*ripz-riz*ripy
1911             gy=riz*ripx-rix*ripz
1912             gz=rix*ripy-riy*ripx
1913             ggm=gmx*gx+gmy*gy+gmz*gz
1914             gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1915             dg=sqrt(gx**2+gy**2+gz**2)
1916             dg3=dg**3
1917             drip=dis(i+1,i+2)
1918             gcr_x=gy*riz-gz*riy
1919             gcr_y=gz*rix-gx*riz
1920             gcr_z=gx*riy-gy*rix
1921             d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1922             d_gcr3=d_gcr**3
1923             gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1924      &           -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1925             if(dgm3.gt.0.0D0.and.
1926      &           dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1927      &           ) then
1928                d20=1.0D0/(dg*dgm)
1929                d21=ggm/(dgm3*dg)
1930                d22=ggm/(dgm*dg3)
1931                d23=1.0D0/(dgm*drip)
1932                d24=gmrp/(dgm3*drip)
1933                sd20=1.0D0/(d_gcr*dgm)
1934                sd21=gcr__gm/(d_gcr3*dgm)
1935                sd22=gcr__gm/(d_gcr*dgm3)
1936             else
1937                d20=0.0D0
1938                d21=0.0D0
1939                d22=0.0D0
1940                d23=0.0D0
1941                d24=0.0D0
1942                sd20=0.0D0
1943                sd21=0.0D0
1944                sd22=0.0D0
1945             endif
1946             atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1947      &           -(ciy*gmz-ciz*gmy)*d21
1948      &           +(ripy*gz-ripz*gy)*d22
1949             atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1950      &           -(ciz*gmx-cix*gmz)*d21
1951      &           +(ripz*gx-ripx*gz)*d22
1952             atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1953      &           -(cix*gmy-ciy*gmx)*d21
1954      &           +(ripx*gy-ripy*gx)*d22
1955 cc**********************************************************************
1956             astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1957      &           -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1958      &           +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1959      &           -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1960      &           +gcr_z*(-ripz*rix+gy))
1961      &           -sd22*(-gmy*ciz+gmz*ciy)
1962
1963             astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1964      &           +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1965      &           +riz*ripz*gmy)
1966      &           -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1967      &           -gcr_z*(ripz*riy+gx))
1968      &           -sd22*(gmx*ciz-gmz*cix)
1969
1970             astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1971      &           +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1972      &           -riz*gx*cix)
1973      &           -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1974      &           +gcr_z*(ripy*riy+ripx*rix))
1975      &           -sd22*(-gmx*ciy+gmy*cix)
1976 cc**********************************************************************
1977 c
1978             apmx(i)=(ciy*ripz-ripy*ciz)*d23
1979      &           -(ciy*gmz-ciz*gmy)*d24
1980             apmy(i)=(ciz*ripx-ripz*cix)*d23
1981      &           -(ciz*gmx-cix*gmz)*d24
1982             apmz(i)=(cix*ripy-ripx*ciy)*d23
1983      &           -(cix*gmy-ciy*gmx)*d24
1984
1985          endif
1986       enddo
1987
1988       return
1989       end
1990 c     END of angvectors
1991 c-------------------------------------------------------------------------------
1992 C---------------------------------------------------------------------------------
1993       subroutine sheetforce(nca,wshet)
1994       implicit none
1995 C     JYLEE
1996 c     this should be matched with dfa.fcm
1997       integer maxca
1998       parameter(maxca=800)
1999 cc**********************************************************************
2000       integer nca
2001       integer i,k
2002       integer inb,nmax,iselect
2003
2004 c      real*8 dfaexp(15001)
2005
2006       real*8 vbeta,vbetp,vbetm
2007       real*8 shefx(maxca,12)
2008       real*8 shefy(maxca,12),shefz(maxca,12)
2009       real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
2010       real*8 vbet(maxca,maxca)
2011       real*8 wshet(maxca,maxca)
2012       real*8 bx(maxca),by(maxca),bz(maxca)
2013
2014       common /sheca/  bx,by,bz
2015       common /phys1/  inb,nmax,iselect
2016       common /shef/   shefx,shefy,shefz
2017       common /shee/   vbeta,vbet,vbetp,vbetm
2018       common /shetf/  shetfx,shetfy,shetfz
2019
2020       inb=nca
2021       do i=1,inb
2022          shetfx(i)=0.0D0
2023          shetfy(i)=0.0D0
2024          shetfz(i)=0.0D0
2025       enddo
2026
2027       do k=1,12
2028          do i=1,inb
2029             shefx(i,k)=0.0D0
2030             shefy(i,k)=0.0D0
2031             shefz(i,k)=0.0D0
2032          enddo
2033       enddo
2034
2035       call sheetene(nca,wshet)
2036       call sheetforce1
2037
2038  887  format(a,1x,i6,3x,f12.8)
2039  888  format(a,1x,i4,1x,i4,3x,f12.8)
2040  889  format(a,1x,i4,3x,f12.8)
2041       !write(2,*) 'coord : '
2042       do i=1,inb
2043          !write(2,887) 'bx:',i,bx(i)
2044          !write(2,887) 'by:',i,by(i)
2045          !write(2,887) 'bz:',i,bz(i)
2046       enddo
2047       !write(2,*) 'After sheetforce1'
2048       do i=1,inb
2049          do k=1,12
2050             !write(2,888) 'shefx :',i,k,shefx(i,k)
2051             !write(2,888) 'shefy :',i,k,shefy(i,k)
2052             !write(2,888) 'shefz :',i,k,shefz(i,k)
2053          enddo
2054       enddo
2055
2056       call sheetforce5
2057
2058       !write(2,*) 'After sheetforce5'
2059       do i=1,inb
2060          do k=1,12
2061             !write(2,888) 'shefx :',i,k,shefx(i,k)
2062             !write(2,888) 'shefy :',i,k,shefy(i,k)
2063             !write(2,888) 'shefz :',i,k,shefz(i,k)
2064          enddo
2065       enddo
2066
2067       call sheetforce6
2068
2069       !write(2,*) 'After sheetforce6'
2070       do i=1,inb
2071          do k=1,12
2072             !write(2,888) 'shefx :',i,k,shefx(i,k)
2073             !write(2,888) 'shefy :',i,k,shefy(i,k)
2074             !write(2,888) 'shefz :',i,k,shefz(i,k)
2075          enddo
2076       enddo
2077
2078       call sheetforce11
2079
2080       !write(2,*) 'After sheetforce11'
2081       do i=1,inb
2082          do k=1,12
2083             !write(2,888) 'shefx :',i,k,shefx(i,k)
2084             !write(2,888) 'shefy :',i,k,shefy(i,k)
2085             !write(2,888) 'shefz :',i,k,shefz(i,k)
2086          enddo
2087       enddo
2088
2089       call sheetforce12
2090
2091       !write(2,*) 'After sheetforce12'
2092       do i=1,inb
2093          do k=1,12
2094             !write(2,888) 'shefx :',i,k,shefx(i,k)
2095             !write(2,888) 'shefy :',i,k,shefy(i,k)
2096             !write(2,888) 'shefz :',i,k,shefz(i,k)
2097          enddo
2098       enddo
2099
2100       do i=1,inb
2101          do k=1,12
2102             shetfx(i)=shetfx(i)+shefx(i,k)
2103             shetfy(i)=shetfy(i)+shefy(i,k)
2104             shetfz(i)=shetfz(i)+shefz(i,k)
2105          enddo
2106       enddo
2107       !write(2,*) 'Beta Finished'
2108       do i=1,inb
2109          !write(2,889) 'shetfx : ',i,shetfx(i)
2110          !write(2,889) 'shetfy : ',i,shetfy(i)
2111          !write(2,889) 'shetfz : ',i,shetfz(i)
2112       enddo
2113
2114       return
2115       end
2116 C     end sheetforce
2117 c-------------------------------------------------------------------------------
2118       subroutine sheetene(nca,wshet)
2119       implicit none
2120       integer maxca
2121       parameter(maxca=800)
2122 cc******************************************************************************
2123
2124 c      real*8 dfaexp(15001)
2125       real*8 dtmp1, dtmp2, dtmp3
2126
2127       real*8 vbet(maxca,maxca)
2128       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2129       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2130       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2131       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2132       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2133       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2134       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2135       real*8 cph(maxca),cth(maxca)
2136       real*8 rx(maxca,maxca)
2137       real*8 ry(maxca,maxca),rz(maxca,maxca)
2138       real*8 bx(maxca),by(maxca),bz(maxca)
2139       real*8 dis(maxca,maxca)
2140       real*8 ulcos(maxca)
2141 cc**********************************************************************
2142       real*8 astx(maxca),asty(maxca),astz(maxca)
2143       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2144       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2145       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2146       real*8 sth(maxca)
2147       real*8 wshet(maxca,maxca)
2148       real*8 dp45, dm45, w_beta
2149       real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2150       integer nca
2151       integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2152       real*8 uum, uup
2153       real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2154
2155       common /sheca/    bx,by,bz
2156       common /phys1/    inb,nmax,iselect
2157       common /kyori2/   dis
2158       common /difvec/   rx,ry,rz
2159       common /coscos/   cph,cth
2160       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2161      &     c00,s00,ulnex,dnex
2162       common /sheconst/ dp45,dm45,w_beta
2163       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2164       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2165       common /shee/    vbeta,vbet,vbetp,vbetm
2166       common /ulang/    ulcos
2167 cc**********************************************************************
2168       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2169      &     astmmz,astm3x,astm3y,astm3z
2170       common /sinsin/   sth
2171
2172       real*8 r_pair_mat(maxca,maxca)
2173 ci      integer istrand(maxca,maxca)
2174 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2175 ci      common  /shetest/ istrand,istrand_p,istrand_m
2176       common /beta_p/ r_pair_mat
2177 C-------------------------------------------------------------------------------
2178       r_pair_mat = 0.0d0
2179       do i=1,inb
2180          do j=1,inb
2181             r_pair_mat(i,j)=wshet(i,j)
2182 c            write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2183          enddo
2184       enddo
2185 c      stop
2186 c
2187       vbeta=0.0D0
2188       vbetp=0.0D0
2189       vbetm=0.0D0
2190
2191       do i=1,inb-7
2192          do j=i+4,inb-3
2193             ip=i+1
2194             ipp=i+2
2195             jp=j+1
2196             jpp=j+2
2197 cc**********************************************************************
2198             y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2199      &           +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2200             y1=-0.5d0*y1/dca
2201             y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2202      &           +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2203             y2=-0.5d0*y2/dnex
2204
2205 cdebug            y2=0
2206
2207             y=y1+y2
2208
2209 ci           if(y.ge.-4) then
2210 ci              istrand(i,j)=1
2211 ci           else
2212 ci              istrand(i,j)=0
2213 ci           endif
2214
2215 ci           if(istrand(i,j).eq.1) then
2216
2217             yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2218             yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2219
2220
2221             pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2222      $           +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2223             pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2224      $           +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2225             pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2226      $           +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2227             pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2228      $           +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2229
2230            yshe1=pin1(i,j)**2+pin2(i,j)**2
2231            yshe1=-0.5d0*yshe1/dshe
2232            yshe2=pin3(i,j)**2+pin4(i,j)**2
2233            yshe2=-0.5d0*yshe2/dshe
2234
2235 ci              if((yshe1+yshe2).ge.-4) then
2236 ci                 istrand_p(i,j)=1
2237 ci              else
2238 ci                 istrand_p(i,j)=0
2239 ci              endif
2240
2241
2242 C            write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2243 C            write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2244 C            write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2245 C            write(*,*) 'dis(i,j):',i,j,dis(i,j)
2246 C            write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2247 C            write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2248 C            write(*,*) 'pin1:',pin1(i,j)
2249 C            write(*,*) 'pin2:',pin2(i,j)
2250 C            write(*,*) 'pin3:',pin3(i,j)
2251 C            write(*,*) 'pin4:',pin4(i,j)
2252
2253 C            write(*,*) 'y:',y
2254 C            write(*,*) 'yy1:',yy1
2255 C            write(*,*) 'yy2:',yy2
2256 C            write(*,*) 'yshe1:',yshe1
2257 C            write(*,*) 'yshe2:',yshe2
2258 c
2259
2260 ci           if (istrand_p(i,j).eq.1) then
2261
2262 cd           yy1=0
2263 cd           yy2=0
2264 cd           yshe1=0
2265 cd           yshe2=0
2266            dtmp1 = y+yy1+yshe1
2267            dtmp2 = y+yy2+yshe2
2268            dtmp3 = y+yy1+yy2+yshe1+yshe2
2269
2270 C            write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2271 C            write(*,*)'2', y,yy1,yy2
2272 C            write(*,*)'3', yshe1,yshe2
2273
2274 cc           if (dtmp3.le.-35.0d0) then
2275 c              vbetap(i,j)=-dp45*exp(dtmp3)
2276 cc              vbetap(i,j)=0.0d0
2277 cc           else
2278 c              vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2279               vbetap(i,j)=-dp45*exp(dtmp3)
2280 cc           end if
2281
2282 cc           if (dtmp1.le.-35.0d0) then
2283 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2284 cc              vbetap1(i,j)=0.0d0
2285 cc           else
2286 c              vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2287 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2288                vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2289 cc           end if
2290
2291 cc           if (dtmp2.le.-35.0d0) then
2292 C              vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2293 cc              vbetap2(i,j)=0.0d0
2294 cc           else
2295 c              vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2296 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2297               vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2298 cc           end if
2299
2300 c           vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2301 c           vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2302 c           vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2303
2304 !           write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2305 !           write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2306
2307 ci           elseif (istrand_p(i,j).eq.0)then
2308 ci            vbetap(i,j)=0
2309 ci            vbetap1(i,j)=0
2310 ci            vbetap2(i,j)=0
2311 ci           endif
2312
2313            yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2314            yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2315
2316            pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2317      $          +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2318            pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2319      $          +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2320            pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2321      $          +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2322            pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2323      $          +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2324
2325            yshe1=pina1(i,j)**2+pina2(i,j)**2
2326            yshe1=-0.5d0*yshe1/dshe
2327            yshe2=pina3(i,j)**2+pina4(i,j)**2
2328            yshe2=-0.5d0*yshe2/dshe
2329
2330 ci              if((yshe1+yshe2).ge.-4) then
2331 ci                 istrand_m(i,j)=1
2332 ci              else
2333 ci                 istrand_m(i,j)=0
2334 ci              endif
2335
2336
2337 C            write(*,*) 'pina1:',pina1(i,j)
2338 C            write(*,*) 'pina2:',pina2(i,j)
2339 C            write(*,*) 'pina3:',pina3(i,j)
2340 C            write(*,*) 'pina4:',pina4(i,j)
2341 C            write(*,*) 'yshe1:',yshe1
2342 C            write(*,*) 'yshe2:',yshe2
2343 C            write(*,*) 'dshe:',dshe
2344
2345 ci           if (istrand_m(i,j).eq.1) then
2346
2347 cd           yy1=0
2348 cd           yy2=0
2349 cd           yshe1=0
2350 cd           yshe2=0
2351
2352            dtmp3=y+yy1+yy2+yshe1+yshe2
2353            dtmp1=y+yy1+yshe1
2354            dtmp2=y+yy2+yshe2
2355
2356 cc           if(dtmp3 .le. -35.0d0) then
2357 c              vbetam(i,j)=-dm45*exp(dtmp3)
2358 cc              vbetam(i,j)=0.0d0
2359 cc           else
2360 c              vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2361               vbetam(i,j)=-dm45*exp(dtmp3)
2362 cc           end if
2363
2364 cc           if(dtmp1 .le. -35.0d0) then
2365 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2366 cc               vbetam1(i,j)=0.0d0
2367 cc           else
2368 c              vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2369 c     $             *dfaexp(idint(-dtmp1*1000)+1)
2370                vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2371 cc           end if
2372
2373 cc           if(dtmp2.le.-35.0d0) then
2374 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2375 cc              vbetam2(i,j)=0.0d0
2376 cc           else
2377 c              vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2378 c     $             *dfaexp(idint(-dtmp2*1000)+1)
2379               vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2380 cc           end if
2381
2382 ci           elseif (istrand_m(i,j).eq.0)then
2383 ci            vbetam(i,j)=0
2384 ci            vbetam1(i,j)=0
2385 ci            vbetam2(i,j)=0
2386 ci           endif
2387
2388
2389 c           vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2390 c           vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2391 c           vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2392
2393 !           write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2394 !           write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2395
2396            uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2397            uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2398
2399 c           write(*,*) 'uup,uum:', uup, uum
2400
2401 c           uup=vbetap1(i,j)+vbetap2(i,j)
2402 c           uum=vbetam1(i,j)+vbetam2(i,j)
2403
2404            vbet(i,j)=uup+uum
2405            vbetp=vbetp+uup
2406            vbetm=vbetm+uum
2407            vbeta=vbeta+vbet(i,j)
2408
2409 ci         elseif(istrand(i,j).eq.0)then
2410 ci           vbet(i,j)=0
2411 ci         endif
2412
2413 c           write(*,*) 'uup,uum:',uup,uum
2414 c           write(*,*) 'vbetap(i,j):',vbetap(i,j)
2415 c           write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2416 c           write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2417 c           write(*,*) 'vbetam(i,j):',vbetam(i,j)
2418 c           write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2419 c           write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2420 c           write(*,*) 'uup:',uup
2421 c           write(*,*) 'uum:',uum
2422 c           write(*,*) 'vbetp:',vbetp
2423 c           write(*,*) 'vbetm:',vbetm
2424 c           write(*,*) 'vbet(i,j):',vbet(i,j)
2425 c           stop
2426
2427         enddo
2428       enddo
2429
2430 !      do i=1,inb-7
2431 !         do j=i+4,inb-3
2432 !            write(*,*) 'I,J:', i,j
2433 !            write(*,*) 'vbetap(i,j):',vbetap(i,j)
2434 !            write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2435 !            write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2436 !            write(*,*) 'vbetam(i,j):',vbetam(i,j)
2437 !            write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2438 !            write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2439 !            write(*,*) 'vbet(i,j):',vbet(i,j)
2440 !         enddo
2441 !      enddo
2442
2443       return
2444       end
2445 c-------------------------------------------------------------------------------
2446       subroutine sheetforce1
2447       implicit none
2448       integer maxca
2449       parameter(maxca=800)
2450 cc**********************************************************************
2451       real*8 vbet(maxca,maxca)
2452       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2453       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2454       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2455       real*8 cph(maxca),cth(maxca)
2456       real*8 rx(maxca,maxca)
2457       real*8 ry(maxca,maxca),rz(maxca,maxca)
2458       real*8 bx(maxca),by(maxca),bz(maxca)
2459       real*8 dis(maxca,maxca)
2460       real*8 shefx(maxca,12)
2461       real*8 shefy(maxca,12),shefz(maxca,12)
2462       real*8 atx(maxca),aty(maxca),atz(maxca)
2463       real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2464       real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2465       real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2466       real*8 apx(maxca),apy(maxca),apz(maxca)
2467       real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2468       real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2469       real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2470       real*8 ulcos(maxca)
2471       real*8 astx(maxca),asty(maxca),astz(maxca)
2472       real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2473       real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2474       real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2475       real*8 sth(maxca)
2476       real*8 w_beta,dp45, dm45
2477       real*8 vbeta, vbetp, vbetm
2478       real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2479      $     c00,s00,ulnex,dnex
2480       integer inb,nmax,iselect
2481
2482       common /phys1/     inb,nmax,iselect
2483       common /kyori2/    dis
2484       common /difvec/   rx,ry,rz
2485       common /coscos/   cph,cth
2486       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2487      $     c00,s00,ulnex,dnex
2488       common /sheconst/ dp45,dm45,w_beta
2489       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2490       common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2491      $     atmmz,atm3x,atm3y,atm3z
2492       common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2493      $     apmmz,apm3x,apm3y,apm3z
2494       common /shef/   shefx,shefy,shefz
2495       common /shee/   vbeta,vbet,vbetp,vbetm
2496       common /ulang/    ulcos
2497 c     c**********************************************************************
2498       common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2499      $     astmmz,astm3x,astm3y,astm3z
2500       common /sinsin/   sth
2501 C--------------------------------------------------------------------------------
2502 c     local variables
2503       integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2504       real*8  c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2505       real*8  c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2506       real*8  c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2507       real*8  dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2508 C--------------------------------------------------------------------------------
2509       do i=4,inb-4
2510          im3=i-3
2511          imm=i-2
2512          im=i-1
2513          c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2514          v1=0.0D0
2515          do j=i+1,inb-3
2516             v1=v1+vbet(im3,j)
2517          enddo
2518          cc1=(ulcos(imm)-ulnex)/dnex
2519          dmm=cc1/(dis(imm,im)*dis(im,i))
2520          dmm__=cc1*ulcos(imm)/dis(im,i)**2
2521          fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2522          fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2523          fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2524 cd         fx=0
2525 cd         fy=0
2526 cd         fz=0
2527          fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2528          fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2529          fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2530          shefx(i,1)=fx*v1
2531          shefy(i,1)=fy*v1
2532          shefz(i,1)=fz*v1
2533       enddo
2534
2535       do i=3,inb-5
2536          imm=i-2
2537          im=i-1
2538          ip=i+1
2539          c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2540          v2=0.0D0
2541          do j=i+2,inb-3
2542             v2=v2+vbet(imm,j)
2543          enddo
2544          cc1=(ulcos(imm)-ulnex)/dnex
2545          cc2=(ulcos(im)-ulnex)/dnex
2546          dmm1=cc1/(dis(imm,im)*dis(im,i))
2547          dmm2=cc2/(dis(im,i)*dis(i,ip))
2548          dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2549          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2550          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2551 cc**********************************************************************
2552          fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2553      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2554          fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2555      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2556          fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2557      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2558 cd         fx=0
2559 cd         fy=0
2560 cd         fz=0
2561          fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2562          fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2563          fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2564          shefx(i,2)=fx*v2
2565          shefy(i,2)=fy*v2
2566          shefz(i,2)=fz*v2
2567       enddo
2568       do i=2,inb-6
2569          im=i-1
2570          ip=i+1
2571          ipp=i+2
2572          c3=(cth(im)*c00+sth(im)*s00-1)/dca
2573          v3=0.0D0
2574          do j=i+3,inb-3
2575             v3=v3+vbet(im,j)
2576          enddo
2577          cc2=(ulcos(im)-ulnex)/dnex
2578          cc3=(ulcos(i)-ulnex)/dnex
2579          dmm2=cc2/(dis(im,i)*dis(i,ip))
2580          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2581          dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2582          dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2583          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2584          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2585      $        -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2586          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2587      $        -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2588          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2589      $        -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2590 cd         fx=0
2591 cd         fy=0
2592 cd         fz=0
2593          fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2594          fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2595          fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2596          shefx(i,3)=fx*v3
2597          shefy(i,3)=fy*v3
2598          shefz(i,3)=fz*v3
2599       enddo
2600       do i=1,inb-7
2601          ip=i+1
2602          ipp=i+2
2603          c4=(cth(i)*c00+sth(i)*s00-1)/dca
2604          v4=0.0D0
2605          do j=i+4,inb-3
2606             v4=v4+vbet(i,j)
2607          enddo
2608          cc3=(ulcos(i)-ulnex)/dnex
2609          dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2610          dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2611          fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2612          fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2613          fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2614 cd         fx=0
2615 cd         fy=0
2616 cd         fz=0
2617          fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2618          fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2619          fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2620          shefx(i,4)=fx*v4
2621          shefy(i,4)=fy*v4
2622          shefz(i,4)=fz*v4
2623       enddo
2624       do j=8,inb
2625          jm3=j-3
2626          jmm=j-2
2627          jm=j-1
2628          c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2629          v7=0.0D0
2630          do i=1,j-7
2631             v7=v7+vbet(i,jm3)
2632          enddo
2633          cc7=(ulcos(jmm)-ulnex)/dnex
2634          dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2635          dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2636          fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2637          fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2638          fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2639 cd         fx=0
2640 cd         fy=0
2641 cd         fz=0
2642          fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2643          fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2644          fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2645          shefx(j,7)=fx*v7
2646          shefy(j,7)=fy*v7
2647          shefz(j,7)=fz*v7
2648       enddo
2649       do j=7,inb-1
2650          jm=j-1
2651          jmm=j-2
2652          jp=j+1
2653          c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2654          v8=0.0D0
2655          do i=1,j-6
2656             v8=v8+vbet(i,jmm)
2657          enddo
2658          cc7=(ulcos(jmm)-ulnex)/dnex
2659          cc8=(ulcos(jm)-ulnex)/dnex
2660          dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2661          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2662          dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2663          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2664          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2665          fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2666      $        -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2667          fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2668      $        -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2669          fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2670      $        -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2671 cd         fx=0
2672 cd         fy=0
2673 cd         fz=0
2674          fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2675          fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2676          fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2677          shefx(j,8)=fx*v8
2678          shefy(j,8)=fy*v8
2679          shefz(j,8)=fz*v8
2680       enddo
2681
2682       do j=6,inb-2
2683          jm=j-1
2684          jp=j+1
2685          jpp=j+2
2686          c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2687          v9=0.0D0
2688          do i=1,j-5
2689             v9=v9+vbet(i,jm)
2690          enddo
2691          cc8=(ulcos(jm)-ulnex)/dnex
2692          cc9=(ulcos(j)-ulnex)/dnex
2693          dmm8=cc8/(dis(jm,j)*dis(j,jp))
2694          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2695          dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2696          dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2697          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2698          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2699      $        -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2700          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2701      $        -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2702          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2703      $        -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2704 cd         fx=0
2705 cd         fy=0
2706 cd         fz=0
2707          fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2708          fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2709          fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2710          shefx(j,9)=fx*v9
2711          shefy(j,9)=fy*v9
2712          shefz(j,9)=fz*v9
2713       enddo
2714
2715       do j=5,inb-3
2716          jp=j+1
2717          jpp=j+2
2718          c10=(cth(j)*c00+sth(j)*s00-1)/dca
2719          v10=0.0D0
2720          do i=1,j-4
2721             v10=v10+vbet(i,j)
2722          enddo
2723          cc9=(ulcos(j)-ulnex)/dnex
2724          dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2725          dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2726          fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2727          fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2728          fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2729 cd         fx=0
2730 cd         fy=0
2731 cd         fz=0
2732          fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2733          fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2734          fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2735          shefx(j,10)=fx*v10
2736          shefy(j,10)=fy*v10
2737          shefz(j,10)=fz*v10
2738       enddo
2739
2740       return
2741       end
2742 c----------------------------------------------------------------------------
2743       subroutine sheetforce5
2744       implicit none
2745       integer maxca
2746       parameter(maxca=800)
2747 cc**********************************************************************
2748       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2749       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2750       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2751       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2752       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2753       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2754       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2755       real*8 rx(maxca,maxca)
2756       real*8 ry(maxca,maxca),rz(maxca,maxca)
2757       real*8 bx(maxca),by(maxca),bz(maxca)
2758       real*8 dis(maxca,maxca)
2759       real*8 shefx(maxca,12),shefy(maxca,12)
2760       real*8 shefz(maxca,12)
2761       real*8 dp45,dm45,w_beta
2762       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2763      $     c00,s00,ulnex,dnex
2764       integer    inb,nmax,iselect
2765 cc**********************************************************************
2766       common /phys1/     inb,nmax,iselect
2767       common /kyori2/    dis
2768       common /difvec/   rx,ry,rz
2769       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2770      $     c00,s00,ulnex,dnex
2771       common /sheconst/ dp45,dm45,w_beta
2772       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2773       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2774       common /shef/   shefx,shefy,shefz
2775 ci      integer istrand(maxca,maxca)
2776 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2777 ci      common  /shetest/ istrand,istrand_p,istrand_m
2778 c********************************************************************************
2779 c     local variables
2780       integer i,imm,im,jp,jpp,j
2781       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2782       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2783       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2784       real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2785       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2786       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2787 c********************************************************************************
2788       do i=3,inb-5
2789          imm=i-2
2790          im=i-1
2791          do j=i+2,inb-3
2792             jp=j+1
2793             jpp=j+2
2794
2795 ci            if(istrand(imm,j).eq.1
2796 ci     &   .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then
2797
2798
2799             yy1=-(dis(i,jpp)-ulhb)/dlhb
2800             y1x=rx(jpp,i)/dis(i,jpp)
2801             y1y=ry(jpp,i)/dis(i,jpp)
2802             y1z=rz(jpp,i)/dis(i,jpp)
2803             y11x=yy1*y1x
2804             y11y=yy1*y1y
2805             y11z=yy1*y1z
2806
2807             yy33=1.0D0/(dis(im,jp)*dis(im,i))
2808             yyy3=pin1(imm,j)/(dis(im,i)**2)
2809             yy3=-pin1(imm,j)/dshe
2810             y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2811             y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2812             y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2813
2814             yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2815             yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2816             yyy4b=pin3(imm,j)/(dis(im,i)**2)
2817             yy4=-pin3(imm,j)/dshe
2818             y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2819      $           -yyy4b*rx(im,i))*yy4
2820             y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2821      $           -yyy4b*ry(im,i))*yy4
2822             y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2823      $           -yyy4b*rz(im,i))*yy4
2824
2825
2826             yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2827             yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2828             yy5=-pin4(imm,j)/dshe
2829             y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2830             y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2831             y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2832
2833             sx=y11x+y3x+y4x+y5x
2834             sy=y11y+y3y+y4y+y5y
2835             sz=y11z+y3z+y4z+y5z
2836
2837             sx1=y3x
2838             sy1=y3y
2839             sz1=y3z
2840             sx2=y11x+y4x+y5x
2841             sy2=y11y+y4y+y5y
2842             sz2=y11z+y4z+y5z
2843
2844             shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2845      $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2846             shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2847      $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2848             shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2849      $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2850
2851 !            shefx(i,5)=shefx(i,5)
2852 !     $           -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2853 !            shefy(i,5)=shefy(i,5)
2854 !     $           -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2855 !            shefz(i,5)=shefz(i,5)
2856 !     $           -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2857
2858             yy6=-(dis(i,jp)-uldhb)/dldhb
2859             y6x=rx(jp,i)/dis(i,jp)
2860             y6y=ry(jp,i)/dis(i,jp)
2861             y6z=rz(jp,i)/dis(i,jp)
2862             y66x=yy6*y6x
2863             y66y=yy6*y6y
2864             y66z=yy6*y6z
2865
2866             yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2867             yyy8=pina1(imm,j)/(dis(im,i)**2)
2868             yy8=-pina1(imm,j)/dshe
2869             y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2870             y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2871             y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2872
2873             yy99=1.0D0/(dis(jp,i)*dis(im,i))
2874             yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2875             yyy9b=pina3(imm,j)/(dis(im,i)**2)
2876             yy9=-pina3(imm,j)/dshe
2877             y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2878      $           -yyy9b*rx(im,i))*yy9
2879             y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2880      $           -yyy9b*ry(im,i))*yy9
2881             y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2882      $           -yyy9b*rz(im,i))*yy9
2883
2884             yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2885             yyy10=pina4(imm,j)/(dis(jp,i)**2)
2886             yy10=-pina4(imm,j)/dshe
2887             y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2888             y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2889             y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2890
2891             sx=y66x+y8x+y9x+y10x
2892             sy=y66y+y8y+y9y+y10y
2893             sz=y66z+y8z+y9z+y10z
2894
2895             sx1=y8x
2896             sy1=y8y
2897             sz1=y8z
2898             sx2=y66x+y9x+y10x
2899             sy2=y66y+y9y+y10y
2900             sz2=y66z+y9z+y10z
2901
2902             shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
2903      $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2904            shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
2905      $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2906             shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
2907      $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2908
2909 !            shefx(i,5)=shefx(i,5)
2910 !     $           -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2911 !            shefy(i,5)=shefy(i,5)
2912 !     $           -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2913 !            shefz(i,5)=shefz(i,5)
2914 !     $           -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2915
2916 ci          endif
2917
2918          enddo
2919       enddo
2920
2921       return
2922       end
2923 c--------------------------------------------------------------------------c
2924       subroutine sheetforce6
2925       implicit none
2926       integer maxca
2927       parameter(maxca=800)
2928 cc**********************************************************************
2929       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2930       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2931       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2932       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2933       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2934       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2935       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2936       real*8 rx(maxca,maxca)
2937       real*8 ry(maxca,maxca),rz(maxca,maxca)
2938       real*8 bx(maxca),by(maxca),bz(maxca)
2939       real*8 dis(maxca,maxca)
2940       real*8 shefx(maxca,12),shefy(maxca,12)
2941       real*8 shefz(maxca,12)
2942       real*8 dp45,dm45,w_beta
2943       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2944      $     c00,s00,ulnex,dnex
2945       integer    inb,nmax,iselect
2946 cc**********************************************************************
2947       common /phys1/     inb,nmax,iselect
2948       common /kyori2/    dis
2949       common /difvec/   rx,ry,rz
2950       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
2951      $     c00,s00,ulnex,dnex
2952       common /sheconst/ dp45,dm45,w_beta
2953       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2954       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2955       common /shef/   shefx,shefy,shefz
2956 ci      integer istrand(maxca,maxca)
2957 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
2958 ci      common  /shetest/ istrand,istrand_p,istrand_m
2959 cc**********************************************************************
2960 C     local variables
2961       integer  i,imm,im,jp,jpp,j,ip
2962       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2963       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2964       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
2965       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2966       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
2967       real*8  yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
2968 C********************************************************************************
2969       do i=2,inb-6
2970          ip=i+1
2971          im=i-1
2972          do j=i+3,inb-3
2973             jp=j+1
2974             jpp=j+2
2975
2976 ci        if(istrand(im,j).eq.1
2977 ci     &    .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then
2978
2979
2980             yy1=-(dis(i,jp)-ulhb)/dlhb
2981             y1x=rx(jp,i)/dis(i,jp)
2982             y1y=ry(jp,i)/dis(i,jp)
2983             y1z=rz(jp,i)/dis(i,jp)
2984             y11x=yy1*y1x
2985             y11y=yy1*y1y
2986             y11z=yy1*y1z
2987
2988             yy33=1.0D0/(dis(i,jp)*dis(i,ip))
2989             yyy3a=pin1(im,j)/(dis(i,jp)**2)
2990             yyy3b=pin1(im,j)/(dis(i,ip)**2)
2991             yy3=-pin1(im,j)/dshe
2992             y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
2993      $           +yyy3b*rx(i,ip))*yy3
2994             y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
2995      $           +yyy3b*ry(i,ip))*yy3
2996             y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
2997      $           +yyy3b*rz(i,ip))*yy3
2998
2999             yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
3000             yyy4=pin2(im,j)/(dis(i,jp)**2)
3001             yy4=-pin2(im,j)/dshe
3002             y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
3003             y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
3004             y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
3005
3006             yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
3007             yyy5=pin3(im,j)/(dis(i,ip)**2)
3008             yy5=-pin3(im,j)/dshe
3009             y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
3010             y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
3011             y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
3012
3013             sx=y11x+y3x+y4x+y5x
3014             sy=y11y+y3y+y4y+y5y
3015             sz=y11z+y3z+y4z+y5z
3016
3017             sx1=y11x+y3x+y4x
3018             sy1=y11y+y3y+y4y
3019             sz1=y11z+y3z+y4z
3020             sx2=y5x
3021             sy2=y5y
3022             sz2=y5z
3023
3024             shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
3025      $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3026             shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
3027      $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3028             shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
3029      $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3030 !            shefx(i,6)=shefx(i,6)
3031 !     $           -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
3032 !            shefy(i,6)=shefy(i,6)
3033 !     $           -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
3034 !            shefz(i,6)=shefz(i,6)
3035 !     $           -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
3036
3037             yy6=-(dis(jpp,i)-uldhb)/dldhb
3038             y6x=rx(jpp,i)/dis(jpp,i)
3039             y6y=ry(jpp,i)/dis(jpp,i)
3040             y6z=rz(jpp,i)/dis(jpp,i)
3041             y66x=yy6*y6x
3042             y66y=yy6*y6y
3043             y66z=yy6*y6z
3044
3045             yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
3046             yyy8a=pina1(im,j)/(dis(i,jpp)**2)
3047             yyy8b=pina1(im,j)/(dis(i,ip)**2)
3048             yy8=-pina1(im,j)/dshe
3049             y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
3050      $           +yyy8b*rx(i,ip))*yy8
3051             y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
3052      $           +yyy8b*ry(i,ip))*yy8
3053             y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
3054      $           +yyy8b*rz(i,ip))*yy8
3055
3056             yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
3057             yyy9=pina2(im,j)/(dis(i,jpp)**2)
3058             yy9=-pina2(im,j)/dshe
3059             y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
3060             y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
3061             y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
3062
3063             yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
3064             yyy10=pina3(im,j)/(dis(i,ip)**2)
3065             yy10=-pina3(im,j)/dshe
3066             y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
3067             y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
3068             y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
3069
3070             sx=y66x+y8x+y9x+y10x
3071             sy=y66y+y8y+y9y+y10y
3072             sz=y66z+y8z+y9z+y10z
3073
3074             sx1=y66x+y8x+y9x
3075             sy1=y66y+y8y+y9y
3076             sz1=y66z+y8z+y9z
3077             sx2=y10x
3078             sy2=y10y
3079             sz2=y10z
3080
3081             shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
3082      $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3083            shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
3084      $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3085             shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
3086      $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3087
3088 !            shefx(i,6)=shefx(i,6)
3089 !     $           -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
3090 !           shefy(i,6)=shefy(i,6)
3091 !     $           -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
3092 !            shefz(i,6)=shefz(i,6)
3093 !     $           -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
3094
3095 ci         endif
3096
3097          enddo
3098       enddo
3099
3100       return
3101       end
3102 c-----------------------------------------------------------------------
3103       subroutine sheetforce11
3104       implicit none
3105       integer maxca
3106       parameter(maxca=800)
3107 cc**********************************************************************
3108       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3109       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3110       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3111       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3112       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3113       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3114       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3115       real*8 rx(maxca,maxca)
3116       real*8 ry(maxca,maxca),rz(maxca,maxca)
3117       real*8 bx(maxca),by(maxca),bz(maxca)
3118       real*8 dis(maxca,maxca)
3119       real*8 shefx(maxca,12),shefy(maxca,12)
3120       real*8 shefz(maxca,12)
3121       real*8 dp45,dm45,w_beta
3122       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3123      $     c00,s00,ulnex,dnex
3124       integer    inb,nmax,iselect
3125 cc**********************************************************************
3126       common /phys1/     inb,nmax,iselect
3127       common /kyori2/    dis
3128       common /difvec/   rx,ry,rz
3129       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3130      $     c00,s00,ulnex,dnex
3131       common /sheconst/ dp45,dm45,w_beta
3132       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3133       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3134       common /shef/   shefx,shefy,shefz
3135 ci      integer istrand(maxca,maxca)
3136 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3137 ci      common  /shetest/ istrand,istrand_p,istrand_m
3138 C********************************************************************************
3139 C     local variables
3140       integer  j,jm,jmm,ip,i,ipp
3141       real*8  yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3142       real*8  yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
3143       real*8  sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
3144       real*8  yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
3145       real*8  yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
3146       real*8  yyy9a,yyy9b,y5z,y66z,y9z,yyy8
3147 C********************************************************************************
3148
3149       do j=7,inb-1
3150          jm=j-1
3151          jmm=j-2
3152          do i=1,j-6
3153             ip=i+1
3154             ipp=i+2
3155
3156 ci            if(istrand(i,jmm).eq.1
3157 ci     &   .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then
3158
3159
3160             yy1=-(dis(ipp,j)-ulhb)/dlhb
3161             y1x=rx(ipp,j)/dis(ipp,j)
3162             y1y=ry(ipp,j)/dis(ipp,j)
3163             y1z=rz(ipp,j)/dis(ipp,j)
3164             y11x=yy1*y1x
3165             y11y=yy1*y1y
3166             y11z=yy1*y1z
3167
3168             yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
3169             yyy3=pin2(i,jmm)/(dis(jm,j)**2)
3170             yy3=-pin2(i,jmm)/dshe
3171             y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
3172             y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
3173             y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
3174
3175             yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
3176             yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
3177             yy4=-pin3(i,jmm)/dshe
3178             y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
3179             y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
3180             y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
3181
3182             yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
3183             yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
3184             yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
3185             yy5=-pin4(i,jmm)/dshe
3186             y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
3187      $           -yyy5b*rx(jm,j))*yy5
3188             y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
3189      $           -yyy5b*ry(jm,j))*yy5
3190             y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
3191      $           -yyy5b*rz(jm,j))*yy5
3192
3193             sx=y11x+y3x+y4x+y5x
3194             sy=y11y+y3y+y4y+y5y
3195             sz=y11z+y3z+y4z+y5z
3196
3197             sx1=y3x
3198             sy1=y3y
3199             sz1=y3z
3200             sx2=y11x+y4x+y5x
3201             sy2=y11y+y4y+y5y
3202             sz2=y11z+y4z+y5z
3203
3204             shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
3205      $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3206             shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
3207      $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3208             shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
3209      $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3210
3211 !            shefx(j,11)=shefx(j,11)
3212 !     $           -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
3213 !            shefy(j,11)=shefy(j,11)
3214 !     $           -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
3215 !            shefz(j,11)=shefz(j,11)
3216 !     $           -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
3217
3218             yy6=-(dis(ip,j)-uldhb)/dldhb
3219             y6x=rx(ip,j)/dis(ip,j)
3220             y6y=ry(ip,j)/dis(ip,j)
3221             y6z=rz(ip,j)/dis(ip,j)
3222             y66x=yy6*y6x
3223             y66y=yy6*y6y
3224             y66z=yy6*y6z
3225
3226             yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
3227             yyy8=pina1(i,jmm)/(dis(ip,j)**2)
3228             yy8=-pina1(i,jmm)/dshe
3229             y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
3230             y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
3231             y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
3232
3233             yy99=1.0D0/(dis(ip,j)*dis(jm,j))
3234             yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
3235             yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
3236             yy9=-pina2(i,jmm)/dshe
3237             y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
3238      $           -yyy9b*rx(jm,j))*yy9
3239             y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
3240      $           -yyy9b*ry(jm,j))*yy9
3241             y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
3242      $           -yyy9b*rz(jm,j))*yy9
3243
3244             yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
3245             yyy10=pina4(i,jmm)/(dis(jm,j)**2)
3246             yy10=-pina4(i,jmm)/dshe
3247             y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
3248             y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
3249             y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
3250
3251             sx=y66x+y8x+y9x+y10x
3252             sy=y66y+y8y+y9y+y10y
3253             sz=y66z+y8z+y9z+y10z
3254
3255             sx1=y66x+y8x+y9x
3256             sy1=y66y+y8y+y9y
3257             sz1=y66z+y8z+y9z
3258             sx2=y10x
3259             sy2=y10y
3260             sz2=y10z
3261
3262             shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3263      $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3264            shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3265      $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3266             shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3267      $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3268
3269 !            shefx(j,11)=shefx(j,11)
3270 !     $           -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3271 !            shefy(j,11)=shefy(j,11)
3272 !     $           -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3273 !            shefz(j,11)=shefz(j,11)
3274 !     $           -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3275
3276 ci         endif
3277
3278          enddo
3279       enddo
3280
3281       return
3282       end
3283 c-----------------------------------------------------------------------
3284       subroutine sheetforce12
3285       implicit none
3286       integer maxca
3287       parameter(maxca=800)
3288 cc**********************************************************************
3289       real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3290       real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3291       real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3292       real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3293       real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3294       real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3295       real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3296       real*8 rx(maxca,maxca)
3297       real*8 ry(maxca,maxca),rz(maxca,maxca)
3298       real*8 bx(maxca),by(maxca),bz(maxca)
3299       real*8 dis(maxca,maxca)
3300       real*8 shefx(maxca,12),shefy(maxca,12)
3301       real*8 shefz(maxca,12)
3302       real*8 dp45,dm45,w_beta
3303       real*8  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3304      $     c00,s00,ulnex,dnex
3305       integer    inb,nmax,iselect
3306 cc**********************************************************************
3307       common /phys1/     inb,nmax,iselect
3308       common /kyori2/    dis
3309       common /difvec/   rx,ry,rz
3310       common /sheparm/  dca,dlhb,ulhb,dshe,dldhb,uldhb,
3311      $     c00,s00,ulnex,dnex
3312       common /sheconst/ dp45,dm45,w_beta
3313       common /she/     vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3314       common /shepin/  pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3315       common /shef/   shefx,shefy,shefz
3316 ci      integer istrand(maxca,maxca)
3317 ci      integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
3318 ci      common  /shetest/ istrand,istrand_p,istrand_m
3319 cc**********************************************************************
3320 C     local variables
3321       integer j,jm,jmm,ip,i,ipp,jp
3322       real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3323       real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3324       real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3325       real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3326       real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3327 !c*************************************************************************c
3328       do j=6,inb-2
3329          jp=j+1
3330          jm=j-1
3331          do i=1,j-5
3332             ip=i+1
3333             ipp=i+2
3334
3335 ci            if(istrand(i,jm).eq.1
3336 ci     &   .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then
3337
3338
3339             yy1=-(dis(ip,j)-ulhb)/dlhb
3340             y1x=rx(ip,j)/dis(ip,j)
3341             y1y=ry(ip,j)/dis(ip,j)
3342             y1z=rz(ip,j)/dis(ip,j)
3343             y11x=y1x*yy1
3344             y11y=y1y*yy1
3345             y11z=y1z*yy1
3346
3347             yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3348             yyy3=pin1(i,jm)/(dis(ip,j)**2)
3349             yy3=-pin1(i,jm)/dshe
3350             y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3351             y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3352             y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3353             yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3354
3355             yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3356             yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3357             yy4=-pin2(i,jm)/dshe
3358             y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3359      $           +yyy4b*rx(j,jp))*yy4
3360             y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3361      $           +yyy4b*ry(j,jp))*yy4
3362             y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3363      $           +yyy4b*rz(j,jp))*yy4
3364
3365             yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3366             yyy5=pin4(i,jm)/(dis(j,jp)**2)
3367             yy5=-pin4(i,jm)/dshe
3368             y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3369             y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3370             y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3371
3372             sx=y11x+y3x+y4x+y5x
3373             sy=y11y+y3y+y4y+y5y
3374             sz=y11z+y3z+y4z+y5z
3375
3376             sx1=y11x+y3x+y4x
3377             sy1=y11y+y3y+y4y
3378             sz1=y11z+y3z+y4z
3379             sx2=y5x
3380             sy2=y5y
3381             sz2=y5z
3382
3383             shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3384      $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3385             shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3386      $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3387             shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3388      $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3389
3390 !            shefx(j,12)=shefx(j,12)
3391 !     $           -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3392 !            shefy(j,12)=shefy(j,12)
3393 !     $           -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3394 !            shefz(j,12)=shefz(j,12)
3395 !     $           -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3396
3397             yy6=-(dis(ipp,j)-uldhb)/dldhb
3398             y6x=rx(ipp,j)/dis(ipp,j)
3399             y6y=ry(ipp,j)/dis(ipp,j)
3400             y6z=rz(ipp,j)/dis(ipp,j)
3401             y66x=yy6*y6x
3402             y66y=yy6*y6y
3403             y66z=yy6*y6z
3404
3405             yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3406             yyy8=pina2(i,jm)/(dis(j,jp)**2)
3407             yy8=-pina2(i,jm)/dshe
3408             y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3409             y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3410             y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3411
3412             yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3413             yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3414             yy9=-pina3(i,jm)/dshe
3415             y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3416             y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3417             y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3418
3419             yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3420             yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3421             yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3422             yy10=-pina4(i,jm)/dshe
3423             y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3424      $           +yyy10b*rx(j,jp))*yy10
3425             y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3426      $           +yyy10b*ry(j,jp))*yy10
3427             y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3428      $           +yyy10b*rz(j,jp))*yy10
3429
3430             sx=y66x+y8x+y9x+y10x
3431             sy=y66y+y8y+y9y+y10y
3432             sz=y66z+y8z+y9z+y10z
3433
3434             sx1=y8x
3435             sy1=y8y
3436             sz1=y8z
3437             sx2=y66x+y9x+y10x
3438             sy2=y66y+y9y+y10y
3439             sz2=y66z+y9z+y10z
3440
3441             shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3442      $           -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3443            shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3444      $           -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3445             shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3446      $           -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3447
3448 ci         endif
3449
3450          ENDDO
3451       ENDDO
3452
3453       RETURN
3454       END
3455 C===============================================================================