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