source/cluster/wham/src/dfa.F

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