source/unres/src_CSA/dfa__.F

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