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