source/unres/src_CSA/dfa.F.not_gly

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