source/unres/src_CSA/dfa.F.org

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