1 subroutine init_dfa_vars
27 CC WEIGHTS for each min
49 C weights for each E term
50 C these should be identical with
57 C precalculate exp table!
60 dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
67 subroutine read_dfa_info
69 C read fragment informations
71 implicit real*8 (a-h,o-z)
73 include 'COMMON.IOUNITS'
74 include 'COMMON.CHAIN'
78 C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
79 C THIS SHOULD BE MODIFIED!!
86 integer ica1, ica2,ica3,ica4,ica5
87 integer ishell, inca, itmp,iitmp
92 open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
94 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
97 write(iout,'(a)') 'dist_dfa.dat is opened!'
99 read(iodfa, '(a)') buffer
100 C read number of restraints
101 read(iodfa, '(i)') IDFADIS
102 read(iodfa, *) dis_inc
104 read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
122 write(iout,'(a)') 'dist_dfa.dat is closed!'
124 C READ ANGLE RESTRAINTS
126 open(iodfa, file='phi_dfa.dat',status='old',err=35)
128 35 write(iout,'(a)') 'Error opening dist_dfa.dat file'
132 write(iout,'(a)') 'phi_dfa.dat is opened!'
135 read(iodfa, '(a)') buffer
136 C READ NUMBER OF RESTRAINTS
137 READ(iodfa, '(i)') IDFAPHI
138 read(iodfa,*) phi_inc
140 read(iodfa,'(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
151 read(iodfa,*) tmp1,tmp2
165 open(iodfa, file='theta_dfa.dat',status='old',err=41)
167 41 write(iout,'(a)') 'Error opening dist_dfa.dat file'
170 write(iout,'(a)') 'theta_dfa.dat is opened!'
172 read(iodfa, '(a)') buffer
173 C READ NUMBER OF RESTRAINTS
174 READ(iodfa, '(i)') IDFATHE
175 read(iodfa,*) the_inc
178 read(iodfa, '(5(i10,1x),1x,i10)')ica1,ica2,ica3,ica4,ica5,nval
189 read(iodfa,*) tmp1,tmp2
201 C END of READING ANGLE RESTRAINT!
203 C NUMBER OF NEIGHBOR CAs
204 open(iodfa,file='nei_dfa.dat',status='old',err=37)
206 37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
209 write(iout,'(a)') 'nei_dfa.dat is opened!'
211 read(iodfa, '(a)') buffer
212 C READ NUMBER OF RESTRAINTS
213 READ(iodfa, '(i)') idfanei
214 read(iodfa,*) nei_inc
217 read(iodfa,'(2(i10,1x),i10)')ica1,ishell,nval
226 C write(*,*) 'READ NEI:',i,j,fnei(i,j)
236 C END OF NEIGHBORING CA
238 C READ BETA RESTRAINT
239 open(iodfa, file='beta_dfa.dat',status='old',err=39)
241 39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
244 write(iout,'(a)') 'beta_dfa.dat is opened!'
246 read(iodfa,'(a)') buffer
247 read(iodfa,'(i)') itmp
248 read(iodfa,*) beta_inc
251 read(iodfa,*) ica1, iitmp
255 c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j)
260 C END OF BETA RESTRAINT
265 subroutine edfad(edfadis)
267 implicit real*8 (a-h,o-z)
269 include 'COMMON.CHAIN'
270 include 'COMMON.DERIV'
273 double precision edfadis
274 integer i, iatm1, iatm2, idiff
275 double precision ckk, sckk,dist,texp
276 double precision jix,jiy,jiz,ep,fp,scc
284 idiff=abs(iatm1-iatm2)
286 jix=c(1,iatm2)-c(1,iatm1)
287 jiy=c(2,iatm2)-c(2,iatm1)
288 jiz=c(3,iatm2)-c(3,iatm1)
289 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
298 write(*,*) 'DIST:', i, dist
299 write(*,*) 'JIX,JIY,JIZ:', jix, jiy, jiz
305 if (dtmp.ge.15.0d0) then
308 texp = dfaexp( idint(dtmp*1000)+1 )/sckk
311 ep=ep+sccdist(i,j)*texp
312 fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
315 write(*,*) 'i,j:', i,j
316 write(*,*) 'iatm1,iatm2:', iatm1, iatm2
317 write(*,*) 'dd,fdist:', dd, fdist(i,j)
318 write(*,*) 'ck,sck:', ckk, sckk
319 write(*,*) 'ep,fp:', ep, fp
320 write(*,*) 'scc:', scc
327 if(abs(ep).lt.1.0d-20)then
333 write(*,*) 'ep,fp:', ep, fp
334 write(*,*) 'edfadis:', edfadis
335 write(*,*) 'gdfad:', gdfad(:,iatm1)
336 write(*,*) 'gdfad:', gdfad(:,iatm2)
338 edfadis = edfadis + ep*dis_inc
340 gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc
341 gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc
342 gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc
344 gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc
345 gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc
346 gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc
348 write(*,*) 'ep,fp:', ep, fp
349 write(*,*) 'edfadis:', edfadis
350 write(*,*) 'gdfad:', gdfad(:,iatm1)
351 write(*,*) 'gdfad:', gdfad(:,iatm2)
353 C gdfad(1,iatm1) = gdfad(1,iatm1)+jix/dist*fp*dis_inc*wwdist
354 C gdfad(2,iatm1) = gdfad(2,iatm1)+jiy/dist*fp*dis_inc*wwdist
355 C gdfad(3,iatm1) = gdfad(3,iatm1)+jiz/dist*fp*dis_inc*wwdist
357 C gdfad(1,iatm2) = gdfad(1,iatm2)-jix/dist*fp*dis_inc*wwdist
358 C gdfad(2,iatm2) = gdfad(2,iatm2)-jiy/dist*fp*dis_inc*wwdist
359 C gdfad(3,iatm2) = gdfad(3,iatm2)-jiz/dist*fp*dis_inc*wwdist
363 print*, 'EDFAD:', edfadis
367 write(*,'(a,i6,1x,3f12.5)')'DFADC:',i,c(1:3,i)
372 write(*,'(a,i6,1x,3f14.7)')'DFADG:', i, gdfad(1:3,i)
378 subroutine edfat(edfator)
380 implicit real*8 (a-h,o-z)
382 include 'COMMON.CHAIN'
383 include 'COMMON.DERIV'
388 double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
389 double precision cwidth, cwidth2
390 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
402 iatom(iii) = iphilis(iii,i)+1
405 C ANGLE VECTOR CALCULTION
406 RIX=C(1,IATOM(2))-C(1,IATOM(1))
407 RIY=C(2,IATOM(2))-C(2,IATOM(1))
408 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
410 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
411 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
412 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
414 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
415 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
416 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
418 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
419 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
420 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
422 GIX=RIY*RIPZ-RIZ*RIPY
423 GIY=RIZ*RIPX-RIX*RIPZ
424 GIZ=RIX*RIPY-RIY*RIPX
426 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
427 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
428 GIPZ=RIPX*RIPPY-RIPY*RIPPX
430 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
431 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
432 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
434 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
435 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
436 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
438 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
439 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
440 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
442 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
443 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
444 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
445 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
447 C END OF ANGLE VECTOR CALCULTION
449 TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
450 APHI(1)=TDOT/(DGI*DRIPP)
451 TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
452 APHI(2)=TDOT/(DGIP*DRIP3)
460 DDPS1=APHI(1)-FPHI1(i,j)
461 DDPS2=APHI(2)-FPHI2(i,j)
463 DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
465 if (dtmp.ge.15.0d0) then
468 ps_tmp = dfaexp(idint(dtmp*1000)+1)
471 ephi=ephi+sccphi(i,j)*ps_tmp
473 tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
474 tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
477 C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
478 C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
481 ephi=-ephi/scc*phi_inc*wwangle
482 tfphi1=tfphi1/scc*phi_inc*wwangle
483 tfphi2=tfphi2/scc*phi_inc*wwangle
485 IF (ABS(EPHI).LT.1d-20) THEN
488 IF (ABS(TFPHI1).LT.1d-20) THEN
491 IF (ABS(TFPHI2).LT.1d-20) THEN
495 C FORCE DIRECTION CALCULATION
500 DM1=1.0d0/(DGI*DRIPP)
502 GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
503 DM2=GIRPP/(DGI**3*DRIPP)
504 DM3=GIRPP/(DGI*DRIPP**3)
506 DM4=1.0d0/(DGIP*DRIP3)
508 GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
509 DM5=GIRP3/(DGIP**3*DRIP3)
510 DM6=GIRP3/(DGIP*DRIP3**3)
512 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
513 & +( GIZ* RIPY- GIY* RIPZ)*DM2
514 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
515 & +( GIX* RIPZ- GIZ* RIPX)*DM2
516 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
517 & +( GIY* RIPX- GIX* RIPY)*DM2
521 C SECOND ATOM BY PHI1
522 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
523 & -(CIPY*GIZ-CIPZ*GIY)*DM2
524 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
525 & -(CIPZ*GIX-CIPX*GIZ)*DM2
526 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
527 & -(CIPX*GIY-CIPY*GIX)*DM2
531 C SECOND ATOM BY PHI2
533 & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
534 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
536 & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
537 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
539 & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
540 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
542 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
543 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
544 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
545 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
546 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
547 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
553 & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
554 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
556 & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
557 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
559 & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
560 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
561 C FOURTH ATOM BY PHI1
562 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
563 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
564 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
565 C FOURTH ATOM BY PHI2
567 & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
568 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
569 & + RIP3X*DM6)*TFPHI2
571 & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
572 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
573 & + RIP3Y*DM6)*TFPHI2
575 & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
576 & -( GIPX*RIPY-RIPX*GIPY)*DM5
577 & + RIP3Z*DM6)*TFPHI2
579 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
580 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
581 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
582 C END OF FORCE DIRECTION
585 gdfat(1,iatom(ii))=gdfat(1,iatom(ii))+tdx(ii)
586 gdfat(2,iatom(ii))=gdfat(2,iatom(ii))+tdy(ii)
587 gdfat(3,iatom(ii))=gdfat(3,iatom(ii))+tdz(ii)
590 enephi = enephi + ephi
591 c end of single assignment statement
593 C END OF PHI RESTRAINT
595 C START OF THETA ANGLE
600 iatom(iii) = iphilis(iii,i)+1
602 c iatom(1:5)=ithelis(1:5,i)
604 C ANGLE VECTOR CALCULTION
605 RIX=C(1,IATOM(2))-C(1,IATOM(1))
606 RIY=C(2,IATOM(2))-C(2,IATOM(1))
607 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
609 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
610 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
611 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
613 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
614 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
615 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
617 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
618 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
619 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
621 GIX=RIY*RIPZ-RIZ*RIPY
622 GIY=RIZ*RIPX-RIX*RIPZ
623 GIZ=RIX*RIPY-RIY*RIPX
625 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
626 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
627 GIPZ=RIPX*RIPPY-RIPY*RIPPX
629 GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
630 GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
631 GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
633 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
634 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
635 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
637 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
638 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
639 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
641 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
642 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
643 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
645 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
646 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
647 DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
648 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
649 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
650 C END OF ANGLE VECTOR CALCULTION
652 TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
653 ATHE(1)=TDOT/(DGI*DGIP)
654 TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
655 ATHE(2)=TDOT/(DGIP*DGIPP)
664 ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
665 ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
666 dtmp= (ddth1**2+ddth2**2)/cwidth2
667 if ( dtmp .ge. 15.0d0) then
670 th_tmp = dfaexp ( idint(dtmp*1000)+1 )
673 ethe=ethe+sccthe(i,j)*th_tmp
675 tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
676 tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
678 C write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
679 C & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
682 ethe=-ethe/scc*the_inc*wwangle
683 tfthe1=tfthe1/scc*the_inc*wwangle
684 tfthe2=tfthe2/scc*the_inc*wwangle
686 if (abs(ethe).lt.tenm20) then
689 if (abs(tfthe1).lt.tenm20) then
692 if (abs(tfthe2).lt.tenm20) then
701 DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
702 DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
704 DM4=1.0d0/(DGIP*DGIPP)
705 DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
706 DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
708 C FIRST ATOM BY THETA1
709 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
710 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
711 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
712 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
713 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
714 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
715 C SECOND ATOM BY THETA1
716 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
717 & -(CIPY*GIZ-CIPZ*GIY)*DM2
718 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
719 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
720 & -(CIPZ*GIX-CIPX*GIZ)*DM2
721 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
722 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
723 & -(CIPX*GIY-CIPY*GIX)*DM2
724 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
725 C SECOND ATOM BY THETA2
727 & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
728 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
730 & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
731 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
733 & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
734 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
735 C THIRD ATOM BY THETA1
736 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
737 & -(GIY*RIZ-GIZ*RIY)*DM2
738 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
739 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
740 & -(GIZ*RIX-GIX*RIZ)*DM2
741 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
742 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
743 & -(GIX*RIY-GIY*RIX)*DM2
744 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
745 C THIRD ATOM BY THETA2
747 & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
748 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
749 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
751 & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
752 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
753 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
755 & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
756 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
757 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
758 C FOURTH ATOM BY THETA1
759 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
760 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
761 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
762 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
763 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
764 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
765 C FOURTH ATOM BY THETA2
767 & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
768 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
769 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
771 & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
772 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
773 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
775 & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
776 & -(GIPX*RIPY-GIPY*RIPX)*DM5
777 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
778 C FIFTH ATOM BY THETA2
779 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
780 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
781 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
782 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
783 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
784 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
785 C !! END OF FORCE DIRECTION!!!!
787 gdfat(1,iatom(ii))=gdfat(1,iatom(ii))+tdx(ii)
788 gdfat(2,iatom(ii))=gdfat(2,iatom(ii))+tdy(ii)
789 gdfat(3,iatom(ii))=gdfat(3,iatom(ii))+tdz(ii)
792 enethe = enethe + ethe
795 edfator = enephi + enethe
797 print*, 'EDFAT:',edfator
800 C write(*,'(a,i6,1x,3f12.5)')'DFATC:',i,c(1:3,i)
805 write(*,'(a,i6,1x,3f14.7)')'DFATG:',i,gdfat(1:3,i)
812 subroutine edfan(edfanei)
813 C DFA neighboring CA restraint
814 implicit real*8 (a-h,o-z)
816 include 'COMMON.CHAIN'
817 include 'COMMON.DERIV'
821 integer kshnum, n1atom
823 double precision enenei,tmp_n
824 double precision pai,hpai
825 double precision jix,jiy,jiz,ndiff,snorm_nei
826 double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
827 double precision dr,dr2,half,ntmp
829 parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
830 parameter(pai=3.14159265358979323846D0)
831 parameter(hpai=1.5707963267948966D0)
832 parameter(snorm_nei=0.886226925452758D0)
838 c print*, 's1:', s1(:)
839 c print*, 's2:', s2(:)
844 n1atom= ineilis(i) + 1
845 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
861 if (n1atom.eq.j) cycle
863 jix=c(1,j)-c(1,n1atom)
864 jiy=c(2,j)-c(2,n1atom)
865 jiz=c(3,j)-c(3,n1atom)
866 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
868 C write(*,*) 'N1ATOM,J,DIST:', n1atom, j, dist
871 if (dist.lt.s1(kshnum).and.
872 & dist.gt.s2(kshnum-1)) then
876 c write(*,*) 'case1:',tmp_n
885 elseif(dist.ge.s1(kshnum).and.
886 & dist.le.s2(kshnum)) then
888 dnei=(dist-s1(kshnum))/dr2*pai
889 tmp_n=tmp_n + half*(1+cos(dnei))
890 c write(*,*) 'case2:',tmp_n
891 ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
901 elseif(dist.ge.s1(kshnum-1).and.
902 & dist.le.s2(kshnum-1)) then
903 dnei=(dist-s1(kshnum-1))/dr2*pai
904 tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
905 c write(*,*) 'case3:',tmp_n
906 ftmp = hpai*sin(dnei)/dr2/dist
918 elseif(kshnum.eq.1) then
920 if(dist.lt.s1(kshnum))then
923 c write(*,*) 'case4:',tmp_n
931 elseif(dist.ge.s1(kshnum).and.
932 & dist.le.s2(kshnum))then
934 dnei=(dist-s1(kshnum))/dr2*pai
935 tmp_n=tmp_n + half*(1+cos(dnei))
936 c write(*,*) 'case5:',tmp_n
937 ftmp = -hpai*sin(dnei)/dr2/dist
958 ndiff = tmp_n-fnei(i,imin)
961 if (dtmp.ge.15.0d0) then
964 ntmp = dfaexp( idint(dtmp*1000) + 1 )
967 enei = enei + sccnei(i,imin)*ntmp
969 & sccnei(i,imin)*ntmp*ndiff*2.0d0
970 scc=scc+sccnei(i,imin)
972 c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
973 c & fnei(i,imin),sccnei(i,imin),enei,scc
976 enei=-enei/scc*snorm_nei*nei_inc*wwnei
977 tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
979 if (abs(enei).lt.1.0d-20)then
982 if (abs(tmp_fnei).lt.1.0d-20) then
992 t2dx(j)=t2dx(j)*tmp_fnei
993 t2dy(j)=t2dy(j)*tmp_fnei
994 t2dz(j)=t2dz(j)*tmp_fnei
997 gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
998 gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
999 gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
1002 gdfan(1,j)=gdfan(1,j)+t2dx(j)
1003 gdfan(2,j)=gdfan(2,j)+t2dy(j)
1004 gdfan(3,j)=gdfan(3,j)+t2dz(j)
1006 c energy calculation
1014 print*, 'NRES: ',nres
1015 print*, 'EDFAN:',edfanei
1018 write(*,'(a,i6,1x,3f12.5)')'DFANC:',i,c(1:3,i)
1023 write(*,'(a,i6,1x,3f14.7)')'DFANG:',i,gdfan(1:3,i)
1029 subroutine edfab(edfabeta)
1031 implicit real*8 (a-h,o-z)
1033 include 'DIMENSIONS'
1034 include 'COMMON.CHAIN'
1035 include 'COMMON.DERIV'
1036 include 'COMMON.DFA'
1039 parameter(pai=3.14159265358979323846d0)
1041 real*8 bx(maxres),by(maxres),bz(maxres)
1042 real*8 vbet(maxres,maxres)
1043 real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
1044 real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
1045 real*8 vbeta,vbetp,vbetm
1046 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
1047 & c00,s00,ulnex,dnex
1048 real*8 dp45,dm45,w_beta
1050 common /sheca/ bx,by,bz
1051 common /shee/ vbeta,vbet,vbetp,vbetm
1052 common /shetf/ shetfx,shetfy,shetfz
1053 common /shef/ shefx, shefy, shefz
1054 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
1055 & c00,s00,ulnex,dnex
1056 common /sheconst/ dp45,dm45,w_beta
1057 C End of sheet variables
1060 double precision enebet
1063 bx=0.0d0;by=0.0d0;bz=0.0d0
1064 shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1079 ulnex=cos(60.0d0/180.0d0*pai)
1086 c00=cos((1.0d0+10.0d0/180.0d0)*pai)
1087 s00=sin((1.0d0+10.0d0/180.0d0)*pai)
1093 C END OF INITIALIZATION
1097 call angvectors(nca)
1098 call sheetforce(nca,wshet,dfaexp)
1100 c end of sheet energy and force
1104 shetfx(j-1)=shetfx(j-1)*beta_inc
1105 shetfy(j-1)=shetfy(j-1)*beta_inc
1106 shetfz(j-1)=shetfz(j-1)*beta_inc
1107 c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1110 vbeta=vbeta*beta_inc
1116 gdfab(1,j)=gdfab(1,j)-shetfx(j-1)
1117 gdfab(2,j)=gdfab(2,j)-shetfy(j-1)
1118 gdfab(3,j)=gdfab(3,j)-shetfz(j-1)
1121 print*, 'EDFAB:',edfabeta
1124 C write(*,'(a,i6,1x,3f12.5)')'DFABC:',i,c(1:3,i)
1127 print*, 'BETA_FORCE'
1129 write(*,'(a,i6,1x,3f14.7)')'DFABG:',i,gdfab(1:3,i)
1134 C-------------------------------------------------------------------------------
1135 subroutine angvectors(nca)
1136 c implicit real*4(a-h,o-z)
1140 parameter(maxca=800)
1142 parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1144 real*8 bx(maxca),by(maxca),bz(maxca)
1145 real*8 dis(maxca,maxca)
1146 real*8 apx(maxca),apy(maxca),apz(maxca)
1147 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
1148 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
1149 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
1150 real*8 atx(maxca),aty(maxca),atz(maxca)
1151 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1152 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1153 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1154 real*8 astx(maxca),asty(maxca),astz(maxca)
1155 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1156 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1157 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1159 real*8 cph(maxca),cth(maxca)
1162 integer i, ip, ipp, ip3, j
1163 real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1164 real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1165 real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1166 real*8 cix, ciy, ciz, cipx, cipy, cipz
1167 real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1168 real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
1169 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1170 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1171 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1172 real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1173 real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1174 real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1175 real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1176 real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1177 real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1178 real*8 grpp,gx,gy,gz
1179 real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1180 real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1181 integer inb,nmax,iselect
1183 common /sheca/ bx,by,bz
1184 common /difvec/ rx, ry, rz
1185 common /ulang/ ulcos
1186 common /phys1/ inb,nmax,iselect
1189 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1190 & apmmz,apm3x,apm3y,apm3z
1191 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1192 & atmmz,atm3x,atm3y,atm3z
1193 common /coscos/ cph,cth
1194 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1195 & astmmz,astm3x,astm3y,astm3z
1197 C-------------------------------------------------------------------------------
1198 c write(*,*) 'inside angvectors'
1203 cph=zero; cth=zero; sth=zero
1204 apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1205 apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1206 atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1207 atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1208 astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1209 astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1212 C r[x,y,z] calc and distance calculation
1213 rx=zero;ry=zero;rz=zero
1220 dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1221 c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1222 c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1223 c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1224 c write(*,*) 'dis(i,j):',i,j,dis(i,j)
1227 c end of r[x,y,z] calc
1233 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1234 ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1235 $ +rz(i,ip)*rz(ip,ipp)
1236 ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1239 c end of virtual bond angle
1240 c write(*,*) 'inside angvectors1'
1251 rippx=bx(ip3)-bx(ipp)
1252 rippy=by(ip3)-by(ipp)
1253 rippz=bz(ip3)-bz(ipp)
1255 gx=riy*ripz-riz*ripy
1256 gy=riz*ripx-rix*ripz
1257 gz=rix*ripy-riy*ripx
1258 gpx=ripy*rippz-ripz*rippy
1259 gpy=ripz*rippx-ripx*rippz
1260 gpz=ripx*rippy-ripy*rippx
1261 gpcrp_x=gpy*ripz-gpz*ripy
1262 gpcrp_y=gpz*ripx-gpx*ripz
1263 gpcrp_z=gpx*ripy-gpy*ripx
1264 d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1265 gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1266 & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1272 gmx=rimy*riz-rimz*riy
1273 gmy=rimz*rix-rimx*riz
1274 gmz=rimx*riy-rimy*rix
1275 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1277 ggm=gmx*gx+gmy*gy+gmz*gz
1278 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1284 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1286 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1287 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1289 c write(*,*) 'inside angvectors2'
1291 rimmx=bx(i-1)-bx(i-2)
1292 rimmy=by(i-1)-by(i-2)
1293 rimmz=bz(i-1)-bz(i-2)
1295 gmmx=rimmy*rimz-rimmz*rimy
1296 gmmy=rimmz*rimx-rimmx*rimz
1297 gmmz=rimmx*rimy-rimmy*rimx
1298 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1300 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1301 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1302 gmcrim_x=gmy*rimz-gmz*rimy
1303 gmcrim_y=gmz*rimx-gmx*rimz
1304 gmcrim_z=gmx*rimy-gmy*rimx
1305 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1306 d_gmcrim3=d_gmcrim**3
1307 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1308 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1312 rim3x=bx(i-2)-bx(i-3)
1313 rim3y=by(i-2)-by(i-3)
1314 rim3z=bz(i-2)-bz(i-3)
1315 g3x=rim3y*rimmz-rim3z*rimmy
1316 g3y=rim3z*rimmx-rim3x*rimmz
1317 g3z=rim3x*rimmy-rim3y*rimmx
1318 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1319 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1320 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1321 cc**********************************************************************
1322 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1323 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1324 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1325 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1326 d_gmmcrimm3=d_gmmcrimm**3
1327 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1328 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1335 dg=sqrt(gx**2+gy**2+gz**2)
1336 dgp=sqrt(gpx**2+gpy**2+gpz**2)
1339 ggp=gx*gpx+gy*gpy+gz*gpz
1340 grpp=gx*rippx+gy*rippy+gz*rippz
1342 if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1343 & .and.d_gpcrp.gt.0.0D0) then
1344 cph(i)=grpp/dg/dripp
1346 sth(i)=gpcrp__g/d_gpcrp/dg
1354 c write(*,*) 'inside angvectors3'
1356 if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1357 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1360 d12=1.0D0/(dg*dripp)
1361 d13=grpp/(dg3*dripp)
1362 sd10=1.0D0/(d_gpcrp*dg)
1363 sd11=gpcrp__g/(d_gpcrp*dg3)
1373 atx(i)=(ripz*gpy-ripy*gpz)*d10
1374 & -(gy*ripz-gz*ripy)*d11
1375 aty(i)=(ripx*gpz-ripz*gpx)*d10
1376 & -(gz*ripx-gx*ripz)*d11
1377 atz(i)=(ripy*gpx-ripx*gpy)*d10
1378 & -(gx*ripy-gy*ripx)*d11
1379 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1380 & +ripy*gpy*ripx-gpx*ripz**2)
1381 & -sd11*(gy*ripz-gz*ripy)
1382 asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1383 & -gpy*ripx**2+gpz*ripy*ripz)
1384 & -sd11*(-gx*ripz+gz*ripx)
1385 astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1386 & -gpz*ripy**2+ripz*gpx*ripx)
1387 & -sd11*(gx*ripy-gy*ripx)
1388 apx(i)=(ripz*rippy-ripy*rippz)*d12
1389 & -(gy*ripz-gz*ripy)*d13
1390 apy(i)=(ripx*rippz-ripz*rippx)*d12
1391 & -(gz*ripx-gx*ripz)*d13
1392 apz(i)=(ripy*rippx-ripx*rippy)*d12
1393 & -(gx*ripy-gy*ripx)*d13
1405 if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1406 & .and.d_gcr3.gt.0.0D0) then
1410 d23=1.0D0/(dgm*drip)
1411 d24=gmrp/(dgm3*drip)
1412 sd20=1.0D0/(d_gcr*dgm)
1413 sd21=gcr__gm/(d_gcr3*dgm)
1414 sd22=gcr__gm/(d_gcr*dgm3)
1425 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1426 & -(ciy*gmz-ciz*gmy)*d21
1427 & +(ripy*gz-ripz*gy)*d22
1428 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1429 & -(ciz*gmx-cix*gmz)*d21
1430 & +(ripz*gx-ripx*gz)*d22
1431 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1432 & -(cix*gmy-ciy*gmx)*d21
1433 & +(ripx*gy-ripy*gx)*d22
1434 cc**********************************************************************
1435 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1436 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1437 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1438 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1439 & +gcr_z*(-ripz*rix+gy))
1440 & -sd22*(-gmy*ciz+gmz*ciy)
1442 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1443 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1445 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1446 & -gcr_z*(ripz*riy+gx))
1447 & -sd22*(gmx*ciz-gmz*cix)
1449 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1450 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1452 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1453 & +gcr_z*(ripy*riy+ripx*rix))
1454 & -sd22*(-gmx*ciy+gmy*cix)
1455 cc**********************************************************************
1456 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1457 & -(ciy*gmz-ciz*gmy)*d24
1458 apmy(i)=(ciz*ripx-ripz*cix)*d23
1459 & -(ciz*gmx-cix*gmz)*d24
1460 apmz(i)=(cix*ripy-ripx*ciy)*d23
1461 & -(cix*gmy-ciy*gmx)*d24
1465 if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1466 & .and.d_gmcrim3.gt.0.0D0) then
1467 d30=1.0D0/(dgm*dgmm)
1468 d31=gmmgm/(dgm3*dgmm)
1469 d32=gmmgm/(dgm*dgmm3)
1470 d33=1.0D0/(dgmm*dri)
1471 d34=gmmr/(dgmm3*dri)
1472 d35=gmmr/(dgmm*dri3)
1473 sd30=1.0D0/(d_gmcrim*dgmm)
1474 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1475 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1488 c write(*,*) 'inside angvectors4'
1490 cc**********************************************************************
1491 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1492 & -(ciy*gmz-ciz*gmy)*d31
1493 & -(gmmy*rimmz-gmmz*rimmy)*d32
1494 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1495 & -(ciz*gmx-cix*gmz)*d31
1496 & -(gmmz*rimmx-gmmx*rimmz)*d32
1497 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1498 & -(cix*gmy-ciy*gmx)*d31
1499 & -(gmmx*rimmy-gmmy*rimmx)*d32
1500 cc**********************************************************************
1501 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1502 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1503 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1504 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1505 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1506 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1507 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1509 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1510 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1511 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1512 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1513 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1514 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1515 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1517 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1518 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1519 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1520 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1521 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1522 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1523 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1524 c**********************************************************************
1525 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1526 & -(gmmy*rimmz-gmmz*rimmy)*d34
1528 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1529 & -(gmmz*rimmx-gmmx*rimmz)*d34
1531 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1532 & -(gmmx*rimmy-gmmy*rimmx)*d34
1537 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1538 & .and.drim3.gt.0.0D0
1539 & .and.d_gmmcrimm3.gt.0.0D0) then
1540 d40=1.0D0/(dg30*dgmm)
1541 d41=g3gmm/(dg30*dgmm3)
1542 d42=1.0D0/(dg30*drim)
1543 d43=g3rim_/(dg30*drim3)
1544 sd40=1.0D0/(dg30*d_gmmcrimm)
1545 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1554 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1555 & -(gmmy*rimmz-gmmz*rimmy)*d41
1556 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1557 & -(gmmz*rimmx-gmmx*rimmz)*d41
1558 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1559 & -(gmmx*rimmy-gmmy*rimmx)*d41
1560 cc**********************************************************************
1561 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1562 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1563 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1564 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1566 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1567 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1568 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1569 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1571 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1572 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1573 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1574 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1575 c**********************************************************************
1576 apm3x(i)=g3x*d42-rimx*d43
1577 apm3y(i)=g3y*d42-rimy*d43
1578 apm3z(i)=g3z*d42-rimz*d43
1581 c*******************************************************************************
1583 c write(*,*) 'inside angvectors5'
1588 rimz=bz(i)-bz(i-1)
\r
1589 rimmx=bx(i-1)-bx(i-2)
1590 rimmy=by(i-1)-by(i-2)
1591 rimmz=bz(i-1)-bz(i-2)
1592 rim3x=bx(i-2)-bx(i-3)
1593 rim3y=by(i-2)-by(i-3)
1594 rim3z=bz(i-2)-bz(i-3)
\r
1595 gmmx=rimmy*rimz-rimmz*rimy
1596 gmmy=rimmz*rimx-rimmx*rimz
1597 gmmz=rimmx*rimy-rimmy*rimx
1598 g3x=rim3y*rimmz-rim3z*rimmy
1599 g3y=rim3z*rimmx-rim3x*rimmz
1600 g3z=rim3x*rimmy-rim3y*rimmx
1602 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1603 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1604 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1609 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1610 cc**********************************************************************
1611 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1612 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1613 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1614 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1615 d_gmmcrimm3=d_gmmcrimm**3
1616 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1617 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1619 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1620 & .and.drim3.gt.0.0D0
1621 & .and.d_gmmcrimm3.gt.0.0D0) then
1622 d40=1.0D0/(dg30*dgmm)
1623 d41=g3gmm/(dg30*dgmm3)
1624 d42=1.0D0/(dg30*drim)
1625 d43=g3rim_/(dg30*drim3)
1626 sd40=1.0D0/(dg30*d_gmmcrimm)
1627 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1636 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1637 & -(gmmy*rimmz-gmmz*rimmy)*d41
1638 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1639 & -(gmmz*rimmx-gmmx*rimmz)*d41
1640 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1641 & -(gmmx*rimmy-gmmy*rimmx)*d41
1642 cc**********************************************************************
1643 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1644 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1645 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1646 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1648 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1649 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1650 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1651 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1653 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1654 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1655 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1656 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1657 cc**********************************************************************
1658 apm3x(i)=g3x*d42-rimx*d43
1659 apm3y(i)=g3y*d42-rimy*d43
1660 apm3z(i)=g3z*d42-rimz*d43
1670 gmx=rimy*riz-rimz*riy
1671 gmy=rimz*rix-rimx*riz
1672 gmz=rimx*riy-rimy*rix
1673 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1677 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1678 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1679 gmcrim_x=gmy*rimz-gmz*rimy
1680 gmcrim_y=gmz*rimx-gmx*rimz
1681 gmcrim_z=gmx*rimy-gmy*rimx
1682 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1683 d_gmcrim3=d_gmcrim**3
1684 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1685 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1687 if(dgm3.gt.0.0D0.and.
1688 & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1689 & .and.d_gmcrim3.gt.0.0D0) then
1690 d30=1.0D0/(dgm*dgmm)
1691 d31=gmmgm/(dgm3*dgmm)
1692 d32=gmmgm/(dgm*dgmm3)
1693 d33=1.0D0/(dgmm*dri)
1694 d34=gmmr/(dgmm3*dri)
1695 d35=gmmr/(dgmm*dri3)
1696 sd30=1.0D0/(d_gmcrim*dgmm)
1697 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1698 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1711 cc**********************************************************************
1712 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1713 & -(ciy*gmz-ciz*gmy)*d31
1714 & -(gmmy*rimmz-gmmz*rimmy)*d32
1715 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1716 & -(ciz*gmx-cix*gmz)*d31
1717 & -(gmmz*rimmx-gmmx*rimmz)*d32
1718 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1719 & -(cix*gmy-ciy*gmx)*d31
1720 & -(gmmx*rimmy-gmmy*rimmx)*d32
1721 cc**********************************************************************
1722 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1723 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1724 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1725 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1726 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1727 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1728 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1730 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1731 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1732 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1733 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1734 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1735 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1736 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1738 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1739 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1740 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1741 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1742 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1743 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1744 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1745 cc**********************************************************************
1746 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1747 & -(gmmy*rimmz-gmmz*rimmy)*d34
1749 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1750 & -(gmmz*rimmx-gmmx*rimmz)*d34
1752 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1753 & -(gmmx*rimmy-gmmy*rimmx)*d34
1757 c write(*,*) 'inside angvectors6'
1767 gx=riy*ripz-riz*ripy
1768 gy=riz*ripx-rix*ripz
1769 gz=rix*ripy-riy*ripx
1770 ggm=gmx*gx+gmy*gy+gmz*gz
1771 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1772 dg=sqrt(gx**2+gy**2+gz**2)
1778 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1780 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1781 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1782 if(dgm3.gt.0.0D0.and.
1783 & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1788 d23=1.0D0/(dgm*drip)
1789 d24=gmrp/(dgm3*drip)
1790 sd20=1.0D0/(d_gcr*dgm)
1791 sd21=gcr__gm/(d_gcr3*dgm)
1792 sd22=gcr__gm/(d_gcr*dgm3)
1803 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1804 & -(ciy*gmz-ciz*gmy)*d21
1805 & +(ripy*gz-ripz*gy)*d22
1806 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1807 & -(ciz*gmx-cix*gmz)*d21
1808 & +(ripz*gx-ripx*gz)*d22
1809 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1810 & -(cix*gmy-ciy*gmx)*d21
1811 & +(ripx*gy-ripy*gx)*d22
1812 cc**********************************************************************
1813 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1814 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1815 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1816 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1817 & +gcr_z*(-ripz*rix+gy))
1818 & -sd22*(-gmy*ciz+gmz*ciy)
1820 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1821 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1823 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1824 & -gcr_z*(ripz*riy+gx))
1825 & -sd22*(gmx*ciz-gmz*cix)
1827 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1828 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1830 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1831 & +gcr_z*(ripy*riy+ripx*rix))
1832 & -sd22*(-gmx*ciy+gmy*cix)
1833 cc**********************************************************************
1835 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1836 & -(ciy*gmz-ciz*gmy)*d24
1837 apmy(i)=(ciz*ripx-ripz*cix)*d23
1838 & -(ciz*gmx-cix*gmz)*d24
1839 apmz(i)=(cix*ripy-ripx*ciy)*d23
1840 & -(cix*gmy-ciy*gmx)*d24
1848 c-------------------------------------------------------------------------------
1849 C---------------------------------------------------------------------------------
1850 subroutine sheetforce(nca,wshet,dfaexp)
1853 c this should be matched with dfa.fcm
1855 parameter(maxca=800)
1856 cc**********************************************************************
1859 integer inb,nmax,iselect
1861 real*8 dfaexp(15001)
1863 real*8 vbeta,vbetp,vbetm
1864 real*8 shefx(maxca,12)
1865 real*8 shefy(maxca,12),shefz(maxca,12)
1866 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
1867 real*8 vbet(maxca,maxca)
1868 real*8 wshet(maxca,maxca)
1869 real*8 bx(maxca),by(maxca),bz(maxca)
1871 common /sheca/ bx,by,bz
1872 common /phys1/ inb,nmax,iselect
1873 common /shef/ shefx,shefy,shefz
1874 common /shee/ vbeta,vbet,vbetp,vbetm
1875 common /shetf/ shetfx,shetfy,shetfz
1892 call sheetene(nca,wshet,dfaexp)
1895 887 format(a,1x,i6,3x,f12.8)
1896 888 format(a,1x,i4,1x,i4,3x,f12.8)
1897 889 format(a,1x,i4,3x,f12.8)
1898 !write(2,*) 'coord : '
1900 !write(2,887) 'bx:',i,bx(i)
1901 !write(2,887) 'by:',i,by(i)
1902 !write(2,887) 'bz:',i,bz(i)
1904 !write(2,*) 'After sheetforce1'
1907 !write(2,888) 'shefx :',i,k,shefx(i,k)
1908 !write(2,888) 'shefy :',i,k,shefy(i,k)
1909 !write(2,888) 'shefz :',i,k,shefz(i,k)
1915 !write(2,*) 'After sheetforce5'
1918 !write(2,888) 'shefx :',i,k,shefx(i,k)
1919 !write(2,888) 'shefy :',i,k,shefy(i,k)
1920 !write(2,888) 'shefz :',i,k,shefz(i,k)
1926 !write(2,*) 'After sheetforce6'
1929 !write(2,888) 'shefx :',i,k,shefx(i,k)
1930 !write(2,888) 'shefy :',i,k,shefy(i,k)
1931 !write(2,888) 'shefz :',i,k,shefz(i,k)
1937 !write(2,*) 'After sheetforce11'
1940 !write(2,888) 'shefx :',i,k,shefx(i,k)
1941 !write(2,888) 'shefy :',i,k,shefy(i,k)
1942 !write(2,888) 'shefz :',i,k,shefz(i,k)
1948 !write(2,*) 'After sheetforce12'
1951 !write(2,888) 'shefx :',i,k,shefx(i,k)
1952 !write(2,888) 'shefy :',i,k,shefy(i,k)
1953 !write(2,888) 'shefz :',i,k,shefz(i,k)
1959 shetfx(i)=shetfx(i)+shefx(i,k)
1960 shetfy(i)=shetfy(i)+shefy(i,k)
1961 shetfz(i)=shetfz(i)+shefz(i,k)
1964 !write(2,*) 'Beta Finished'
1966 !write(2,889) 'shetfx : ',i,shetfx(i)
1967 !write(2,889) 'shetfy : ',i,shetfy(i)
1968 !write(2,889) 'shetfz : ',i,shetfz(i)
1974 c-------------------------------------------------------------------------------
1975 subroutine sheetene(nca,wshet,dfaexp)
1978 parameter(maxca=800)
1979 cc******************************************************************************
1981 real*8 dfaexp(15001)
1982 real*8 dtmp1, dtmp2, dtmp3
1984 real*8 vbet(maxca,maxca)
1985 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
1986 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
1987 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
1988 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
1989 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
1990 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
1991 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
1992 real*8 cph(maxca),cth(maxca)
1993 real*8 rx(maxca,maxca)
1994 real*8 ry(maxca,maxca),rz(maxca,maxca)
1995 real*8 bx(maxca),by(maxca),bz(maxca)
1996 real*8 dis(maxca,maxca)
1998 cc**********************************************************************
1999 real*8 astx(maxca),asty(maxca),astz(maxca)
2000 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2001 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2002 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2004 real*8 wshet(maxca,maxca)
2005 real*8 dp45, dm45, w_beta
2006 real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
2008 integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
2010 real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
2012 common /sheca/ bx,by,bz
2013 common /phys1/ inb,nmax,iselect
2015 common /difvec/ rx,ry,rz
2016 common /coscos/ cph,cth
2017 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2018 & c00,s00,ulnex,dnex
2019 common /sheconst/ dp45,dm45,w_beta
2020 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2021 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2022 common /shee/ vbeta,vbet,vbetp,vbetm
2023 common /ulang/ ulcos
2024 cc**********************************************************************
2025 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2026 & astmmz,astm3x,astm3y,astm3z
2029 real*8 r_pair_mat(maxca,maxca)
2030 common /beta_p/ r_pair_mat
2031 C-------------------------------------------------------------------------------
2035 r_pair_mat(i,j)=wshet(i,j)
2036 c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
2051 cc**********************************************************************
2052 y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
2053 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
2055 y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
2056 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
2060 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
2061 yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
2063 pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
2064 $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
2065 pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
2066 $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
2067 pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
2068 $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
2069 pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
2070 $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2072 yshe1=pin1(i,j)**2+pin2(i,j)**2
2073 yshe1=-0.5d0*yshe1/dshe
2074 yshe2=pin3(i,j)**2+pin4(i,j)**2
2075 yshe2=-0.5d0*yshe2/dshe
2077 C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2078 C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2079 C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2080 C write(*,*) 'dis(i,j):',i,j,dis(i,j)
2081 C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2082 C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2083 C write(*,*) 'pin1:',pin1(i,j)
2084 C write(*,*) 'pin2:',pin2(i,j)
2085 C write(*,*) 'pin3:',pin3(i,j)
2086 C write(*,*) 'pin4:',pin4(i,j)
2089 C write(*,*) 'yy1:',yy1
2090 C write(*,*) 'yy2:',yy2
2091 C write(*,*) 'yshe1:',yshe1
2092 C write(*,*) 'yshe2:',yshe2
2097 dtmp3 = y+yy1+yy2+yshe1+yshe2
2099 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2100 C write(*,*)'2', y,yy1,yy2
2101 C write(*,*)'3', yshe1,yshe2
2103 if (dtmp3.le.-15.0d0) then
2104 c vbetap(i,j)=-dp45*exp(dtmp3)
2107 vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2110 if (dtmp1.le.-15.0d0) then
2111 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2114 vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2115 $ *dfaexp(idint(-dtmp1*1000)+1)
2118 if (dtmp2.le.-15.0d0) then
2119 C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2122 vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2123 $ *dfaexp(idint(-dtmp2*1000)+1)
2126 c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2127 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2128 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2130 ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2131 ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2134 yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2135 yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2137 pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2138 $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2139 pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2140 $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2141 pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2142 $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2143 pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2144 $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2146 yshe1=pina1(i,j)**2+pina2(i,j)**2
2147 yshe1=-0.5d0*yshe1/dshe
2148 yshe2=pina3(i,j)**2+pina4(i,j)**2
2149 yshe2=-0.5d0*yshe2/dshe
2151 C write(*,*) 'pina1:',pina1(i,j)
2152 C write(*,*) 'pina2:',pina2(i,j)
2153 C write(*,*) 'pina3:',pina3(i,j)
2154 C write(*,*) 'pina4:',pina4(i,j)
2155 C write(*,*) 'yshe1:',yshe1
2156 C write(*,*) 'yshe2:',yshe2
2157 C write(*,*) 'dshe:',dshe
2159 dtmp3=y+yy1+yy2+yshe1+yshe2
2163 if(dtmp3 .le. -15.0d0) then
2164 c vbetam(i,j)=-dm45*exp(dtmp3)
2167 vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2170 if(dtmp1 .le. -15.0d0) then
2171 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2174 vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2175 $ *dfaexp(idint(-dtmp1*1000)+1)
2178 if(dtmp2.le.-15.0d0) then
2179 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2182 vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2183 $ *dfaexp(idint(-dtmp2*1000)+1)
2186 c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2187 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2188 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2190 ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2191 ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2193 uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2194 uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2196 c write(*,*) 'uup,uum:', uup, uum
2198 c uup=vbetap1(i,j)+vbetap2(i,j)
2199 c uum=vbetam1(i,j)+vbetam2(i,j)
2204 vbeta=vbeta+vbet(i,j)
2206 c write(*,*) 'uup,uum:',uup,uum
2207 c write(*,*) 'vbetap(i,j):',vbetap(i,j)
2208 c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2209 c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2210 c write(*,*) 'vbetam(i,j):',vbetam(i,j)
2211 c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2212 c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2213 c write(*,*) 'uup:',uup
2214 c write(*,*) 'uum:',uum
2215 c write(*,*) 'vbetp:',vbetp
2216 c write(*,*) 'vbetm:',vbetm
2217 c write(*,*) 'vbet(i,j):',vbet(i,j)
2225 ! write(*,*) 'I,J:', i,j
2226 ! write(*,*) 'vbetap(i,j):',vbetap(i,j)
2227 ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2228 ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2229 ! write(*,*) 'vbetam(i,j):',vbetam(i,j)
2230 ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2231 ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2232 ! write(*,*) 'vbet(i,j):',vbet(i,j)
2238 c-------------------------------------------------------------------------------
2239 subroutine sheetforce1
2242 parameter(maxca=800)
2243 cc**********************************************************************
2244 real*8 vbet(maxca,maxca)
2245 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2246 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2247 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2248 real*8 cph(maxca),cth(maxca)
2249 real*8 rx(maxca,maxca)
2250 real*8 ry(maxca,maxca),rz(maxca,maxca)
2251 real*8 bx(maxca),by(maxca),bz(maxca)
2252 real*8 dis(maxca,maxca)
2253 real*8 shefx(maxca,12)
2254 real*8 shefy(maxca,12),shefz(maxca,12)
2255 real*8 atx(maxca),aty(maxca),atz(maxca)
2256 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2257 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2258 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2259 real*8 apx(maxca),apy(maxca),apz(maxca)
2260 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2261 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2262 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2264 real*8 astx(maxca),asty(maxca),astz(maxca)
2265 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2266 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2267 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2269 real*8 w_beta,dp45, dm45
2270 real*8 vbeta, vbetp, vbetm
2271 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2272 $ c00,s00,ulnex,dnex
2273 integer inb,nmax,iselect
2275 common /phys1/ inb,nmax,iselect
2277 common /difvec/ rx,ry,rz
2278 common /coscos/ cph,cth
2279 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2280 $ c00,s00,ulnex,dnex
2281 common /sheconst/ dp45,dm45,w_beta
2282 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2283 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2284 $ atmmz,atm3x,atm3y,atm3z
2285 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2286 $ apmmz,apm3x,apm3y,apm3z
2287 common /shef/ shefx,shefy,shefz
2288 common /shee/ vbeta,vbet,vbetp,vbetm
2289 common /ulang/ ulcos
2290 c c**********************************************************************
2291 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2292 $ astmmz,astm3x,astm3y,astm3z
2294 C--------------------------------------------------------------------------------
2296 integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2297 real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2298 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2299 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2300 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2301 C--------------------------------------------------------------------------------
2306 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2311 cc1=(ulcos(imm)-ulnex)/dnex
2312 dmm=cc1/(dis(imm,im)*dis(im,i))
2313 dmm__=cc1*ulcos(imm)/dis(im,i)**2
2314 fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2315 fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2316 fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2317 fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2318 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2319 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2329 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2334 cc1=(ulcos(imm)-ulnex)/dnex
2335 cc2=(ulcos(im)-ulnex)/dnex
2336 dmm1=cc1/(dis(imm,im)*dis(im,i))
2337 dmm2=cc2/(dis(im,i)*dis(i,ip))
2338 dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2339 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2340 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2341 cc**********************************************************************
2342 fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2343 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2344 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2345 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2346 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2347 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2348 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2349 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2350 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2359 c3=(cth(im)*c00+sth(im)*s00-1)/dca
2364 cc2=(ulcos(im)-ulnex)/dnex
2365 cc3=(ulcos(i)-ulnex)/dnex
2366 dmm2=cc2/(dis(im,i)*dis(i,ip))
2367 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2368 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2369 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2370 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2371 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2372 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2373 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2374 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2375 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2376 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2377 fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2378 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2379 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2387 c4=(cth(i)*c00+sth(i)*s00-1)/dca
2392 cc3=(ulcos(i)-ulnex)/dnex
2393 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2394 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2395 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2396 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2397 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2398 fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2399 fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2400 fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2409 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2414 cc7=(ulcos(jmm)-ulnex)/dnex
2415 dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2416 dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2417 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2418 fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2419 fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2420 fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2421 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2422 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2431 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2436 cc7=(ulcos(jmm)-ulnex)/dnex
2437 cc8=(ulcos(jm)-ulnex)/dnex
2438 dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2439 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2440 dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2441 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2442 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2443 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2444 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2445 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2446 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2447 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2448 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2449 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2450 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2451 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2461 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2466 cc8=(ulcos(jm)-ulnex)/dnex
2467 cc9=(ulcos(j)-ulnex)/dnex
2468 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2469 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2470 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2471 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2472 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2473 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2474 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2475 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2476 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2477 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2478 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2479 fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2480 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2481 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2490 c10=(cth(j)*c00+sth(j)*s00-1)/dca
2495 cc9=(ulcos(j)-ulnex)/dnex
2496 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2497 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2498 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2499 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2500 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2501 fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2502 fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2503 fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2511 c----------------------------------------------------------------------------
2512 subroutine sheetforce5
2515 parameter(maxca=800)
2516 cc**********************************************************************
2517 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2518 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2519 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2520 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2521 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2522 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2523 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2524 real*8 rx(maxca,maxca)
2525 real*8 ry(maxca,maxca),rz(maxca,maxca)
2526 real*8 bx(maxca),by(maxca),bz(maxca)
2527 real*8 dis(maxca,maxca)
2528 real*8 shefx(maxca,12),shefy(maxca,12)
2529 real*8 shefz(maxca,12)
2530 real*8 dp45,dm45,w_beta
2531 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2532 $ c00,s00,ulnex,dnex
2533 integer inb,nmax,iselect
2534 cc**********************************************************************
2535 common /phys1/ inb,nmax,iselect
2537 common /difvec/ rx,ry,rz
2538 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2539 $ c00,s00,ulnex,dnex
2540 common /sheconst/ dp45,dm45,w_beta
2541 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2542 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2543 common /shef/ shefx,shefy,shefz
2544 c********************************************************************************
2546 integer i,imm,im,jp,jpp,j
2547 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2548 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2549 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2550 real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2551 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2552 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2553 c********************************************************************************
2561 yy1=-(dis(i,jpp)-ulhb)/dlhb
2562 y1x=rx(jpp,i)/dis(i,jpp)
2563 y1y=ry(jpp,i)/dis(i,jpp)
2564 y1z=rz(jpp,i)/dis(i,jpp)
2569 yy33=1.0D0/(dis(im,jp)*dis(im,i))
2570 yyy3=pin1(imm,j)/(dis(im,i)**2)
2571 yy3=-pin1(imm,j)/dshe
2572 y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2573 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2574 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2576 yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2577 yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2578 yyy4b=pin3(imm,j)/(dis(im,i)**2)
2579 yy4=-pin3(imm,j)/dshe
2580 y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2581 $ -yyy4b*rx(im,i))*yy4
2582 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2583 $ -yyy4b*ry(im,i))*yy4
2584 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2585 $ -yyy4b*rz(im,i))*yy4
2588 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2589 yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2590 yy5=-pin4(imm,j)/dshe
2591 y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2592 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2593 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2606 shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2607 $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2608 shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2609 $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2610 shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2611 $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2613 ! shefx(i,5)=shefx(i,5)
2614 ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2615 ! shefy(i,5)=shefy(i,5)
2616 ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2617 ! shefz(i,5)=shefz(i,5)
2618 ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2620 yy6=-(dis(i,jp)-uldhb)/dldhb
2621 y6x=rx(jp,i)/dis(i,jp)
2622 y6y=ry(jp,i)/dis(i,jp)
2623 y6z=rz(jp,i)/dis(i,jp)
2628 yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2629 yyy8=pina1(imm,j)/(dis(im,i)**2)
2630 yy8=-pina1(imm,j)/dshe
2631 y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2632 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2633 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2635 yy99=1.0D0/(dis(jp,i)*dis(im,i))
2636 yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2637 yyy9b=pina3(imm,j)/(dis(im,i)**2)
2638 yy9=-pina3(imm,j)/dshe
2639 y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2640 $ -yyy9b*rx(im,i))*yy9
2641 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2642 $ -yyy9b*ry(im,i))*yy9
2643 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2644 $ -yyy9b*rz(im,i))*yy9
2646 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2647 yyy10=pina4(imm,j)/(dis(jp,i)**2)
2648 yy10=-pina4(imm,j)/dshe
2649 y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2650 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2651 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2653 sx=y66x+y8x+y9x+y10x
2654 sy=y66y+y8y+y9y+y10y
2655 sz=y66z+y8z+y9z+y10z
2664 shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
2665 $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2666 shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
2667 $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2668 shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
2669 $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2671 ! shefx(i,5)=shefx(i,5)
2672 ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2673 ! shefy(i,5)=shefy(i,5)
2674 ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2675 ! shefz(i,5)=shefz(i,5)
2676 ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2683 c--------------------------------------------------------------------------c
2684 subroutine sheetforce6
2687 parameter(maxca=800)
2688 cc**********************************************************************
2689 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2690 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2691 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2692 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2693 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2694 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2695 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2696 real*8 rx(maxca,maxca)
2697 real*8 ry(maxca,maxca),rz(maxca,maxca)
2698 real*8 bx(maxca),by(maxca),bz(maxca)
2699 real*8 dis(maxca,maxca)
2700 real*8 shefx(maxca,12),shefy(maxca,12)
2701 real*8 shefz(maxca,12)
2702 real*8 dp45,dm45,w_beta
2703 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2704 $ c00,s00,ulnex,dnex
2705 integer inb,nmax,iselect
2706 cc**********************************************************************
2707 common /phys1/ inb,nmax,iselect
2709 common /difvec/ rx,ry,rz
2710 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2711 $ c00,s00,ulnex,dnex
2712 common /sheconst/ dp45,dm45,w_beta
2713 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2714 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2715 common /shef/ shefx,shefy,shefz
2716 cc**********************************************************************
2718 integer i,imm,im,jp,jpp,j,ip
2719 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2720 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2721 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
2722 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2723 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
2724 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
2725 C********************************************************************************
2733 yy1=-(dis(i,jp)-ulhb)/dlhb
2734 y1x=rx(jp,i)/dis(i,jp)
2735 y1y=ry(jp,i)/dis(i,jp)
2736 y1z=rz(jp,i)/dis(i,jp)
2741 yy33=1.0D0/(dis(i,jp)*dis(i,ip))
2742 yyy3a=pin1(im,j)/(dis(i,jp)**2)
2743 yyy3b=pin1(im,j)/(dis(i,ip)**2)
2744 yy3=-pin1(im,j)/dshe
2745 y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
2746 $ +yyy3b*rx(i,ip))*yy3
2747 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
2748 $ +yyy3b*ry(i,ip))*yy3
2749 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
2750 $ +yyy3b*rz(i,ip))*yy3
2752 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
2753 yyy4=pin2(im,j)/(dis(i,jp)**2)
2754 yy4=-pin2(im,j)/dshe
2755 y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
2756 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
2757 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
2759 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
2760 yyy5=pin3(im,j)/(dis(i,ip)**2)
2761 yy5=-pin3(im,j)/dshe
2762 y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
2763 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
2764 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
2777 shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
2778 $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
2779 shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
2780 $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
2781 shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
2782 $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
2783 ! shefx(i,6)=shefx(i,6)
2784 ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
2785 ! shefy(i,6)=shefy(i,6)
2786 ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
2787 ! shefz(i,6)=shefz(i,6)
2788 ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
2790 yy6=-(dis(jpp,i)-uldhb)/dldhb
2791 y6x=rx(jpp,i)/dis(jpp,i)
2792 y6y=ry(jpp,i)/dis(jpp,i)
2793 y6z=rz(jpp,i)/dis(jpp,i)
2798 yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
2799 yyy8a=pina1(im,j)/(dis(i,jpp)**2)
2800 yyy8b=pina1(im,j)/(dis(i,ip)**2)
2801 yy8=-pina1(im,j)/dshe
2802 y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
2803 $ +yyy8b*rx(i,ip))*yy8
2804 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
2805 $ +yyy8b*ry(i,ip))*yy8
2806 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
2807 $ +yyy8b*rz(i,ip))*yy8
2809 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2810 yyy9=pina2(im,j)/(dis(i,jpp)**2)
2811 yy9=-pina2(im,j)/dshe
2812 y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
2813 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
2814 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
2816 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
2817 yyy10=pina3(im,j)/(dis(i,ip)**2)
2818 yy10=-pina3(im,j)/dshe
2819 y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
2820 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
2821 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
2823 sx=y66x+y8x+y9x+y10x
2824 sy=y66y+y8y+y9y+y10y
2825 sz=y66z+y8z+y9z+y10z
2834 shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
2835 $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
2836 shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
2837 $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
2838 shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
2839 $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
2841 ! shefx(i,6)=shefx(i,6)
2842 ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
2843 ! shefy(i,6)=shefy(i,6)
2844 ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
2845 ! shefz(i,6)=shefz(i,6)
2846 ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
2853 c-----------------------------------------------------------------------
2854 subroutine sheetforce11
2857 parameter(maxca=800)
2858 cc**********************************************************************
2859 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2860 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2861 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2862 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2863 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2864 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2865 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2866 real*8 rx(maxca,maxca)
2867 real*8 ry(maxca,maxca),rz(maxca,maxca)
2868 real*8 bx(maxca),by(maxca),bz(maxca)
2869 real*8 dis(maxca,maxca)
2870 real*8 shefx(maxca,12),shefy(maxca,12)
2871 real*8 shefz(maxca,12)
2872 real*8 dp45,dm45,w_beta
2873 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2874 $ c00,s00,ulnex,dnex
2875 integer inb,nmax,iselect
2876 cc**********************************************************************
2877 common /phys1/ inb,nmax,iselect
2879 common /difvec/ rx,ry,rz
2880 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2881 $ c00,s00,ulnex,dnex
2882 common /sheconst/ dp45,dm45,w_beta
2883 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2884 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2885 common /shef/ shefx,shefy,shefz
2886 C********************************************************************************
2888 integer j,jm,jmm,ip,i,ipp
2889 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2890 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
2891 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
2892 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
2893 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
2894 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
2895 C********************************************************************************
2904 yy1=-(dis(ipp,j)-ulhb)/dlhb
2905 y1x=rx(ipp,j)/dis(ipp,j)
2906 y1y=ry(ipp,j)/dis(ipp,j)
2907 y1z=rz(ipp,j)/dis(ipp,j)
2912 yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
2913 yyy3=pin2(i,jmm)/(dis(jm,j)**2)
2914 yy3=-pin2(i,jmm)/dshe
2915 y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
2916 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
2917 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
2919 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
2920 yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
2921 yy4=-pin3(i,jmm)/dshe
2922 y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
2923 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
2924 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
2926 yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
2927 yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
2928 yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
2929 yy5=-pin4(i,jmm)/dshe
2930 y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
2931 $ -yyy5b*rx(jm,j))*yy5
2932 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
2933 $ -yyy5b*ry(jm,j))*yy5
2934 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
2935 $ -yyy5b*rz(jm,j))*yy5
2948 shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
2949 $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
2950 shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
2951 $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
2952 shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
2953 $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
2955 ! shefx(j,11)=shefx(j,11)
2956 ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
2957 ! shefy(j,11)=shefy(j,11)
2958 ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
2959 ! shefz(j,11)=shefz(j,11)
2960 ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
2962 yy6=-(dis(ip,j)-uldhb)/dldhb
2963 y6x=rx(ip,j)/dis(ip,j)
2964 y6y=ry(ip,j)/dis(ip,j)
2965 y6z=rz(ip,j)/dis(ip,j)
2970 yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
2971 yyy8=pina1(i,jmm)/(dis(ip,j)**2)
2972 yy8=-pina1(i,jmm)/dshe
2973 y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
2974 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
2975 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
2977 yy99=1.0D0/(dis(ip,j)*dis(jm,j))
2978 yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
2979 yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
2980 yy9=-pina2(i,jmm)/dshe
2981 y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
2982 $ -yyy9b*rx(jm,j))*yy9
2983 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
2984 $ -yyy9b*ry(jm,j))*yy9
2985 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
2986 $ -yyy9b*rz(jm,j))*yy9
2988 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
2989 yyy10=pina4(i,jmm)/(dis(jm,j)**2)
2990 yy10=-pina4(i,jmm)/dshe
2991 y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
2992 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
2993 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
2995 sx=y66x+y8x+y9x+y10x
2996 sy=y66y+y8y+y9y+y10y
2997 sz=y66z+y8z+y9z+y10z
3006 shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
3007 $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3008 shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
3009 $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3010 shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
3011 $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3013 ! shefx(j,11)=shefx(j,11)
3014 ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
3015 ! shefy(j,11)=shefy(j,11)
3016 ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
3017 ! shefz(j,11)=shefz(j,11)
3018 ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
3025 c-----------------------------------------------------------------------
3026 subroutine sheetforce12
3029 parameter(maxca=800)
3030 cc**********************************************************************
3031 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
3032 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
3033 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
3034 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
3035 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
3036 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
3037 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
3038 real*8 rx(maxca,maxca)
3039 real*8 ry(maxca,maxca),rz(maxca,maxca)
3040 real*8 bx(maxca),by(maxca),bz(maxca)
3041 real*8 dis(maxca,maxca)
3042 real*8 shefx(maxca,12),shefy(maxca,12)
3043 real*8 shefz(maxca,12)
3044 real*8 dp45,dm45,w_beta
3045 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
3046 $ c00,s00,ulnex,dnex
3047 integer inb,nmax,iselect
3048 cc**********************************************************************
3049 common /phys1/ inb,nmax,iselect
3051 common /difvec/ rx,ry,rz
3052 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
3053 $ c00,s00,ulnex,dnex
3054 common /sheconst/ dp45,dm45,w_beta
3055 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
3056 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
3057 common /shef/ shefx,shefy,shefz
3058 cc**********************************************************************
3060 integer j,jm,jmm,ip,i,ipp,jp
3061 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
3062 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
3063 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
3064 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3065 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3066 !c*************************************************************************c
3074 yy1=-(dis(ip,j)-ulhb)/dlhb
3075 y1x=rx(ip,j)/dis(ip,j)
3076 y1y=ry(ip,j)/dis(ip,j)
3077 y1z=rz(ip,j)/dis(ip,j)
3082 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3083 yyy3=pin1(i,jm)/(dis(ip,j)**2)
3084 yy3=-pin1(i,jm)/dshe
3085 y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3086 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3087 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3088 yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3090 yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3091 yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3092 yy4=-pin2(i,jm)/dshe
3093 y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3094 $ +yyy4b*rx(j,jp))*yy4
3095 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3096 $ +yyy4b*ry(j,jp))*yy4
3097 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3098 $ +yyy4b*rz(j,jp))*yy4
3100 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3101 yyy5=pin4(i,jm)/(dis(j,jp)**2)
3102 yy5=-pin4(i,jm)/dshe
3103 y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3104 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3105 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3118 shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3119 $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3120 shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3121 $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3122 shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3123 $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3125 ! shefx(j,12)=shefx(j,12)
3126 ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3127 ! shefy(j,12)=shefy(j,12)
3128 ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3129 ! shefz(j,12)=shefz(j,12)
3130 ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3132 yy6=-(dis(ipp,j)-uldhb)/dldhb
3133 y6x=rx(ipp,j)/dis(ipp,j)
3134 y6y=ry(ipp,j)/dis(ipp,j)
3135 y6z=rz(ipp,j)/dis(ipp,j)
3140 yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3141 yyy8=pina2(i,jm)/(dis(j,jp)**2)
3142 yy8=-pina2(i,jm)/dshe
3143 y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3144 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3145 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3147 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3148 yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3149 yy9=-pina3(i,jm)/dshe
3150 y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3151 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3152 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3154 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3155 yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3156 yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3157 yy10=-pina4(i,jm)/dshe
3158 y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3159 $ +yyy10b*rx(j,jp))*yy10
3160 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3161 $ +yyy10b*ry(j,jp))*yy10
3162 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3163 $ +yyy10b*rz(j,jp))*yy10
3165 sx=y66x+y8x+y9x+y10x
3166 sy=y66y+y8y+y9y+y10y
3167 sz=y66z+y8z+y9z+y10z
3176 shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3177 $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3178 shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3179 $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3180 shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3181 $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3188 C===============================================================================