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