source/wham/src-HCD-5D/dfa.F

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