1 subroutine init_dfa_vars
4 include 'COMMON.INTERACT'
28 CC WEIGHTS for each min
45 C weights for each E term
46 C these should be identical with
53 C precalculate exp table!
56 c dfaexp(ii) = exp(-ii*0.001d0 + 0.0005d0)
62 print *,'ishiftca=',ishiftca,'ilastca=',ilastca
68 subroutine read_dfa_info
70 C read fragment informations
72 implicit real*8 (a-h,o-z)
74 include 'COMMON.IOUNITS'
75 include 'COMMON.CHAIN'
79 C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
80 C THIS SHOULD BE MODIFIED!!
87 integer ica1, ica2,ica3,ica4,ica5
88 integer ishell, inca, itmp,iitmp
93 open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
95 33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
98 write(iout,'(a)') 'dist_dfa.dat is opened!'
100 read(iodfa, '(a)') buffer
101 C read number of restraints
102 read(iodfa, '(i)') IDFADIS
103 read(iodfa, *) dis_inc
105 read(iodfa, '(i10,1x,i10,1x,i10)') ica1, ica2, nval
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
282 iatm1=idislis(1,i)+ishiftca
283 iatm2=idislis(2,i)+ishiftca
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)
302 if (dtmp.ge.15.0d0) then
305 c texp = dfaexp( idint(dtmp*1000)+1 )/sckk
306 texp = exp(-dtmp)/sckk
309 ep=ep+sccdist(i,j)*texp
310 fp=fp+sccdist(i,j)*texp*dd*2.0d0/ckk
312 C write(*,'(2i8,6f12.5)') i, j, dist,
313 C & fdist(i,j), ep, fp, sccdist(i,j), scc
321 c IF(ABS(EP).lt.1.0d-20)THEN
324 c IF (ABS(FP).lt.1.0d-20) THEN
328 edfadis=edfadis+ep*dis_inc*wwdist
330 gdfad(1,iatm1) = gdfad(1,iatm1)-jix/dist*fp*dis_inc*wwdist
331 gdfad(2,iatm1) = gdfad(2,iatm1)-jiy/dist*fp*dis_inc*wwdist
332 gdfad(3,iatm1) = gdfad(3,iatm1)-jiz/dist*fp*dis_inc*wwdist
334 gdfad(1,iatm2) = gdfad(1,iatm2)+jix/dist*fp*dis_inc*wwdist
335 gdfad(2,iatm2) = gdfad(2,iatm2)+jiy/dist*fp*dis_inc*wwdist
336 gdfad(3,iatm2) = gdfad(3,iatm2)+jiz/dist*fp*dis_inc*wwdist
343 subroutine edfat(edfator)
345 implicit real*8 (a-h,o-z)
347 include 'COMMON.CHAIN'
348 include 'COMMON.DERIV'
353 double precision aphi(2),athe(2),tdx(5),tdy(5),tdz(5)
354 double precision cwidth, cwidth2
355 PARAMETER(CWIDTH=0.1D0,CWIDTH2=0.2D0,PAI=3.14159265358979323846D0)
367 iatom(iii)=iphilis(iii,i)+ishiftca
370 C ANGLE VECTOR CALCULTION
371 RIX=C(1,IATOM(2))-C(1,IATOM(1))
372 RIY=C(2,IATOM(2))-C(2,IATOM(1))
373 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
375 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
376 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
377 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
379 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
380 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
381 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
383 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
384 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
385 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
387 GIX=RIY*RIPZ-RIZ*RIPY
388 GIY=RIZ*RIPX-RIX*RIPZ
389 GIZ=RIX*RIPY-RIY*RIPX
391 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
392 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
393 GIPZ=RIPX*RIPPY-RIPY*RIPPX
395 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
396 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
397 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
399 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
400 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
401 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
403 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
404 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
405 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
407 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
408 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
409 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
410 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
412 C END OF ANGLE VECTOR CALCULTION
414 TDOT=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
415 APHI(1)=TDOT/(DGI*DRIPP)
416 TDOT=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
417 APHI(2)=TDOT/(DGIP*DRIP3)
425 DDPS1=APHI(1)-FPHI1(i,j)
426 DDPS2=APHI(2)-FPHI2(i,j)
428 DTMP = (DDPS1**2+DDPS2**2)/CWIDTH2
430 if (dtmp.ge.15.0d0) then
433 c ps_tmp = dfaexp(idint(dtmp*1000)+1)
437 ephi=ephi+sccphi(i,j)*ps_tmp
439 tfphi1=tfphi1+sccphi(i,j)*ddps1/cwidth*ps_tmp
440 tfphi2=tfphi2+sccphi(i,j)*ddps2/cwidth*ps_tmp
443 C write(*,'(2i8,8f12.6)')i,j,aphi(1),fphi1(i,j),
444 C & aphi(2),fphi2(i,j),tfphi1,tfphi2,ephi,sccphi(i,j)
447 ephi=-ephi/scc*phi_inc*wwangle
448 tfphi1=tfphi1/scc*phi_inc*wwangle
449 tfphi2=tfphi2/scc*phi_inc*wwangle
451 IF (ABS(EPHI).LT.1d-20) THEN
454 IF (ABS(TFPHI1).LT.1d-20) THEN
457 IF (ABS(TFPHI2).LT.1d-20) THEN
461 C FORCE DIRECTION CALCULATION
466 DM1=1.0d0/(DGI*DRIPP)
468 GIRPP=GIX*RIPPX+GIY*RIPPY+GIZ*RIPPZ
469 DM2=GIRPP/(DGI**3*DRIPP)
470 DM3=GIRPP/(DGI*DRIPP**3)
472 DM4=1.0d0/(DGIP*DRIP3)
474 GIRP3=GIPX*RIP3X+GIPY*RIP3Y+GIPZ*RIP3Z
475 DM5=GIRP3/(DGIP**3*DRIP3)
476 DM6=GIRP3/(DGIP*DRIP3**3)
478 TDX(1)=(RIPZ*RIPPY-RIPY*RIPPZ)*DM1
479 & +( GIZ* RIPY- GIY* RIPZ)*DM2
480 TDY(1)=(RIPX*RIPPZ-RIPZ*RIPPX)*DM1
481 & +( GIX* RIPZ- GIZ* RIPX)*DM2
482 TDZ(1)=(RIPY*RIPPX-RIPX*RIPPY)*DM1
483 & +( GIY* RIPX- GIX* RIPY)*DM2
487 C SECOND ATOM BY PHI1
488 TDX(2)=(CIPY*RIPPZ-CIPZ*RIPPY)*DM1
489 & -(CIPY*GIZ-CIPZ*GIY)*DM2
490 TDY(2)=(CIPZ*RIPPX-CIPX*RIPPZ)*DM1
491 & -(CIPZ*GIX-CIPX*GIZ)*DM2
492 TDZ(2)=(CIPX*RIPPY-CIPY*RIPPX)*DM1
493 & -(CIPX*GIY-CIPY*GIX)*DM2
497 C SECOND ATOM BY PHI2
499 & ((RIPPZ*RIP3Y-RIPPY*RIP3Z)*DM4
500 & +( GIPZ*RIPPY- GIPY*RIPPZ)*DM5)*TFPHI2
502 & ((RIPPX*RIP3Z-RIPPZ*RIP3X)*DM4
503 & +( GIPX*RIPPZ- GIPZ*RIPPX)*DM5)*TFPHI2
505 & ((RIPPY*RIP3X-RIPPX*RIP3Y)*DM4
506 & +( GIPY*RIPPX- GIPX*RIPPY)*DM5)*TFPHI2
508 TDX(3)=(-GIX+RIPPY*RIZ-RIPPZ*RIY)*DM1
509 & -(GIY*RIZ-RIY*GIZ)*DM2+RIPPX*DM3
510 TDY(3)=(-GIY+RIPPZ*RIX-RIPPX*RIZ)*DM1
511 & -(GIZ*RIX-RIZ*GIX)*DM2+RIPPY*DM3
512 TDZ(3)=(-GIZ+RIPPX*RIY-RIPPY*RIX)*DM1
513 & -(GIX*RIY-RIX*GIY)*DM2+RIPPZ*DM3
519 & ((CIPPY*RIP3Z-CIPPZ*RIP3Y)*DM4
520 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5)*TFPHI2
522 & ((CIPPZ*RIP3X-CIPPX*RIP3Z)*DM4
523 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5)*TFPHI2
525 & ((CIPPX*RIP3Y-CIPPY*RIP3X)*DM4
526 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5)*TFPHI2
527 C FOURTH ATOM BY PHI1
528 TDX(4)=(GIX*DM1-RIPPX*DM3)*TFPHI1
529 TDY(4)=(GIY*DM1-RIPPY*DM3)*TFPHI1
530 TDZ(4)=(GIZ*DM1-RIPPZ*DM3)*TFPHI1
531 C FOURTH ATOM BY PHI2
533 & ((-GIPX+RIP3Y*RIPZ-RIP3Z*RIPY)*DM4
534 & -( GIPY*RIPZ-RIPY*GIPZ)*DM5
535 & + RIP3X*DM6)*TFPHI2
537 & ((-GIPY+RIP3Z*RIPX-RIP3X*RIPZ)*DM4
538 & -( GIPZ*RIPX-RIPZ*GIPX)*DM5
539 & + RIP3Y*DM6)*TFPHI2
541 & ((-GIPZ+RIP3X*RIPY-RIP3Y*RIPX)*DM4
542 & -( GIPX*RIPY-RIPX*GIPY)*DM5
543 & + RIP3Z*DM6)*TFPHI2
545 TDX(5)=(GIPX*DM4-RIP3X*DM6)*TFPHI2
546 TDY(5)=(GIPY*DM4-RIP3Y*DM6)*TFPHI2
547 TDZ(5)=(GIPZ*DM4-RIP3Z*DM6)*TFPHI2
548 C END OF FORCE DIRECTION
551 gdfat(1,IATOM(II))=gdfat(1,IATOM(II))+TDX(II)
552 gdfat(2,IATOM(II))=gdfat(2,IATOM(II))+TDY(II)
553 gdfat(3,IATOM(II))=gdfat(3,IATOM(II))+TDZ(II)
556 enephi = enephi + ephi
557 c end of single assignment statement
559 C END OF PHI RESTRAINT
561 C START OF THETA ANGLE
566 iatom(iii)=ithelis(iii,i)+ishiftca
570 C ANGLE VECTOR CALCULTION
571 RIX=C(1,IATOM(2))-C(1,IATOM(1))
572 RIY=C(2,IATOM(2))-C(2,IATOM(1))
573 RIZ=C(3,IATOM(2))-C(3,IATOM(1))
575 RIPX=C(1,IATOM(3))-C(1,IATOM(2))
576 RIPY=C(2,IATOM(3))-C(2,IATOM(2))
577 RIPZ=C(3,IATOM(3))-C(3,IATOM(2))
579 RIPPX=C(1,IATOM(4))-C(1,IATOM(3))
580 RIPPY=C(2,IATOM(4))-C(2,IATOM(3))
581 RIPPZ=C(3,IATOM(4))-C(3,IATOM(3))
583 RIP3X=C(1,IATOM(5))-C(1,IATOM(4))
584 RIP3Y=C(2,IATOM(5))-C(2,IATOM(4))
585 RIP3Z=C(3,IATOM(5))-C(3,IATOM(4))
587 GIX=RIY*RIPZ-RIZ*RIPY
588 GIY=RIZ*RIPX-RIX*RIPZ
589 GIZ=RIX*RIPY-RIY*RIPX
591 GIPX=RIPY*RIPPZ-RIPZ*RIPPY
592 GIPY=RIPZ*RIPPX-RIPX*RIPPZ
593 GIPZ=RIPX*RIPPY-RIPY*RIPPX
595 GIPPX=RIPPY*RIP3Z-RIPPZ*RIP3Y
596 GIPPY=RIPPZ*RIP3X-RIPPX*RIP3Z
597 GIPPZ=RIPPX*RIP3Y-RIPPY*RIP3X
599 CIPX=C(1,IATOM(3))-C(1,IATOM(1))
600 CIPY=C(2,IATOM(3))-C(2,IATOM(1))
601 CIPZ=C(3,IATOM(3))-C(3,IATOM(1))
603 CIPPX=C(1,IATOM(4))-C(1,IATOM(2))
604 CIPPY=C(2,IATOM(4))-C(2,IATOM(2))
605 CIPPZ=C(3,IATOM(4))-C(3,IATOM(2))
607 CIP3X=C(1,IATOM(5))-C(1,IATOM(3))
608 CIP3Y=C(2,IATOM(5))-C(2,IATOM(3))
609 CIP3Z=C(3,IATOM(5))-C(3,IATOM(3))
611 DGI=SQRT(GIX*GIX+GIY*GIY+GIZ*GIZ)
612 DGIP=SQRT(GIPX*GIPX+GIPY*GIPY+GIPZ*GIPZ)
613 DGIPP=SQRT(GIPPX*GIPPX+GIPPY*GIPPY+GIPPZ*GIPPZ)
614 DRIPP=SQRT(RIPPX*RIPPX+RIPPY*RIPPY+RIPPZ*RIPPZ)
615 DRIP3=SQRT(RIP3X*RIP3X+RIP3Y*RIP3Y+RIP3Z*RIP3Z)
616 C END OF ANGLE VECTOR CALCULTION
618 TDOT=GIX*GIPX+GIY*GIPY+GIZ*GIPZ
619 ATHE(1)=TDOT/(DGI*DGIP)
620 TDOT=GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ
621 ATHE(2)=TDOT/(DGIP*DGIPP)
630 ddth1=athe(1)-fthe1(i,j) !cos(the1)-cos(the1_ref)
631 ddth2=athe(2)-fthe2(i,j) !cos(the2)-cos(the2_ref)
632 dtmp= (ddth1**2+ddth2**2)/cwidth2
633 if ( dtmp .ge. 15.0d0) then
636 c th_tmp = dfaexp ( idint(dtmp*1000)+1 )
640 ethe=ethe+sccthe(i,j)*th_tmp
642 tfthe1=tfthe1+sccthe(i,j)*ddth1/cwidth*th_tmp !-dv/dcos(the1)
643 tfthe2=tfthe2+sccthe(i,j)*ddth2/cwidth*th_tmp !-dv/dcos(the2)
645 C write(*,'(2i8,8f12.6)')i,j,athe(1),fthe1(i,j),
646 C & athe(2),fthe2(i,j),tfthe1,tfthe2,ethe,sccthe(i,j)
649 ethe=-ethe/scc*the_inc*wwangle
650 tfthe1=tfthe1/scc*the_inc*wwangle
651 tfthe2=tfthe2/scc*the_inc*wwangle
653 IF (ABS(ETHE).LT.TENM20) THEN
656 IF (ABS(TFTHE1).LT.TENM20) THEN
659 IF (ABS(TFTHE2).LT.TENM20) THEN
668 DM2=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI**3*DGIP)
669 DM3=(GIX*GIPX+GIY*GIPY+GIZ*GIPZ)/(DGI*DGIP**3)
671 DM4=1.0d0/(DGIP*DGIPP)
672 DM5=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP**3*DGIPP)
673 DM6=(GIPX*GIPPX+GIPY*GIPPY+GIPZ*GIPPZ)/(DGIP*DGIPP**3)
675 C FIRST ATOM BY THETA1
676 TDX(1)=((RIPZ*GIPY-RIPY*GIPZ)*DM1
677 & -(GIY*RIPZ-GIZ*RIPY)*DM2)*TFTHE1
678 TDY(1)=((-RIPZ*GIPX+RIPX*GIPZ)*DM1
679 & -(-GIX*RIPZ+GIZ*RIPX)*DM2)*TFTHE1
680 TDZ(1)=((RIPY*GIPX-RIPX*GIPY)*DM1
681 & -(GIX*RIPY-GIY*RIPX)*DM2)*TFTHE1
682 C SECOND ATOM BY THETA1
683 TDX(2)=((CIPY*GIPZ-CIPZ*GIPY-RIPPY*GIZ+RIPPZ*GIY)*DM1
684 & -(CIPY*GIZ-CIPZ*GIY)*DM2
685 & +(RIPPY*GIPZ-RIPPZ*GIPY)*DM3)*TFTHE1
686 TDY(2)=((CIPZ*GIPX-CIPX*GIPZ-RIPPZ*GIX+RIPPX*GIZ)*DM1
687 & -(CIPZ*GIX-CIPX*GIZ)*DM2
688 & +(RIPPZ*GIPX-RIPPX*GIPZ)*DM3)*TFTHE1
689 TDZ(2)=((CIPX*GIPY-CIPY*GIPX-RIPPX*GIY+RIPPY*GIX)*DM1
690 & -(CIPX*GIY-CIPY*GIX)*DM2
691 & +(RIPPX*GIPY-RIPPY*GIPX)*DM3)*TFTHE1
692 C SECOND ATOM BY THETA2
694 & ((RIPPZ*GIPPY-RIPPY*GIPPZ)*DM4
695 & -(GIPY*RIPPZ-GIPZ*RIPPY)*DM5)*TFTHE2
697 & ((-RIPPZ*GIPPX+RIPPX*GIPPZ)*DM4
698 & -(-GIPX*RIPPZ+GIPZ*RIPPX)*DM5)*TFTHE2
700 & ((RIPPY*GIPPX-RIPPX*GIPPY)*DM4
701 & -(GIPX*RIPPY-GIPY*RIPPX)*DM5)*TFTHE2
702 C THIRD ATOM BY THETA1
703 TDX(3)=((GIPY*RIZ-GIPZ*RIY-GIY*CIPPZ+GIZ*CIPPY)*DM1
704 & -(GIY*RIZ-GIZ*RIY)*DM2
705 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM3) *TFTHE1
706 TDY(3)=((GIPZ*RIX-GIPX*RIZ-GIZ*CIPPX+GIX*CIPPZ)*DM1
707 & -(GIZ*RIX-GIX*RIZ)*DM2
708 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM3) *TFTHE1
709 TDZ(3)=((GIPX*RIY-GIPY*RIX-GIX*CIPPY+GIY*CIPPX)*DM1
710 & -(GIX*RIY-GIY*RIX)*DM2
711 & -(CIPPX*GIPY-CIPPY*GIPX)*DM3) *TFTHE1
712 C THIRD ATOM BY THETA2
714 & ((CIPPY*GIPPZ-CIPPZ*GIPPY-RIP3Y*GIPZ+RIP3Z*GIPY)*DM4
715 & -(CIPPY*GIPZ-CIPPZ*GIPY)*DM5
716 & +(RIP3Y*GIPpZ-RIP3Z*GIPpY)*DM6) *TFTHE2
718 & ((CIPPZ*GIPPX-CIPPX*GIPPZ-RIP3Z*GIPX+RIP3X*GIPZ)*DM4
719 & -(CIPPZ*GIPX-CIPPX*GIPZ)*DM5
720 & +(RIP3Z*GIPpX-RIP3X*GIPpZ)*DM6) *TFTHE2
722 & ((CIPPX*GIPPY-CIPPY*GIPPX-RIP3X*GIPY+RIP3Y*GIPX)*DM4
723 & -(CIPPX*GIPY-CIPPY*GIPX)*DM5
724 & +(RIP3X*GIPpY-RIP3Y*GIPpX)*DM6) *TFTHE2
725 C FOURTH ATOM BY THETA1
726 TDX(4)=-((GIZ*RIPY-GIY*RIPZ)*DM1
727 & -(GIPZ*RIPY-GIPY*RIPZ)*DM3) *TFTHE1
728 TDY(4)=-((GIX*RIPZ-GIZ*RIPX)*DM1
729 & -(GIPX*RIPZ-GIPZ*RIPX)*DM3) *TFTHE1
730 TDZ(4)=-((GIY*RIPX-GIX*RIPY)*DM1
731 & -(GIPY*RIPX-GIPX*RIPY)*DM3) *TFTHE1
732 C FOURTH ATOM BY THETA2
734 & ((GIPPY*RIPZ-GIPPZ*RIPY-GIPY*CIP3Z+GIPZ*CIP3Y)*DM4
735 & -(GIPY*RIPZ-GIPZ*RIPY)*DM5
736 & -(CIP3Y*GIPPZ-CIP3Z*GIPPY)*DM6)*TFTHE2
738 & ((GIPPZ*RIPX-GIPPX*RIPZ-GIPZ*CIP3X+GIPX*CIP3Z)*DM4
739 & -(GIPZ*RIPX-GIPX*RIPZ)*DM5
740 & -(CIP3Z*GIPPX-CIP3X*GIPPZ)*DM6)*TFTHE2
742 & ((GIPPX*RIPY-GIPPY*RIPX-GIPX*CIP3Y+GIPY*CIP3X)*DM4
743 & -(GIPX*RIPY-GIPY*RIPX)*DM5
744 & -(CIP3X*GIPPY-CIP3Y*GIPPX)*DM6)*TFTHE2
745 C FIFTH ATOM BY THETA2
746 TDX(5)=-((GIPZ*RIPPY-GIPY*RIPPZ)*DM4
747 & -(GIPPZ*RIPPY-GIPPY*RIPPZ)*DM6)*TFTHE2
748 TDY(5)=-((GIPX*RIPPZ-GIPZ*RIPPX)*DM4
749 & -(GIPPX*RIPPZ-GIPPZ*RIPPX)*DM6)*TFTHE2
750 TDZ(5)=-((GIPY*RIPPX-GIPX*RIPPY)*DM4
751 & -(GIPPY*RIPPX-GIPPX*RIPPY)*DM6)*TFTHE2
752 C !! END OF FORCE DIRECTION!!!!
754 gdfat(1,iatom(II))=gdfat(1,iatom(II))+TDX(II)
755 gdfat(2,iatom(II))=gdfat(2,iatom(II))+TDY(II)
756 gdfat(3,iatom(II))=gdfat(3,iatom(II))+TDZ(II)
759 enethe = enethe + ethe
762 edfator = enephi + enethe
767 subroutine edfan(edfanei)
768 C DFA neighboring CA restraint
769 implicit real*8 (a-h,o-z)
771 include 'COMMON.CHAIN'
772 include 'COMMON.DERIV'
776 integer kshnum, n1atom
778 double precision enenei,tmp_n
779 double precision pai,hpai
780 double precision jix,jiy,jiz,ndiff,snorm_nei
781 double precision t2dx(maxres),t2dy(maxres),t2dz(maxres)
782 double precision dr,dr2,half,ntmp,dtmp
784 parameter(dr=0.25d0,dr2=0.50d0,half=0.50d0)
785 parameter(pai=3.14159265358979323846D0)
786 parameter(hpai=1.5707963267948966D0)
787 parameter(snorm_nei=0.886226925452758D0)
793 c print*, 's1:', s1(:)
794 c print*, 's2:', s2(:)
799 n1atom=ineilis(i)+ishiftca
800 C write(*,*) 'kshnum,n1atom:', kshnum, n1atom
813 do j = ishiftca+1, ilastca
815 if (n1atom.eq.j) cycle
817 jix=c(1,j)-c(1,n1atom)
818 jiy=c(2,j)-c(2,n1atom)
819 jiz=c(3,j)-c(3,n1atom)
820 dist=sqrt(jix*jix+jiy*jiy+jiz*jiz)
822 c write(*,*) n1atom, j, dist
825 if (dist.lt.s1(kshnum).and.
826 & dist.gt.s2(kshnum-1)) then
830 c write(*,*) 'case1:',tmp_n
839 elseif(dist.ge.s1(kshnum).and.
840 & dist.le.s2(kshnum)) then
842 dnei=(dist-s1(kshnum))/dr2*pai
843 tmp_n=tmp_n + half*(1+cos(dnei))
844 c write(*,*) 'case2:',tmp_n
845 ftmp=-pai*sin(dnei)/dr2/dist/2.0d0
855 elseif(dist.ge.s1(kshnum-1).and.
856 & dist.le.s2(kshnum-1)) then
857 dnei=(dist-s1(kshnum-1))/dr2*pai
858 tmp_n=tmp_n + 1.0d0 - half*(1+cos(dnei))
859 c write(*,*) 'case3:',tmp_n
860 ftmp = hpai*sin(dnei)/dr2/dist
872 elseif(kshnum.eq.1) then
874 if(dist.lt.s1(kshnum))then
877 c write(*,*) 'case4:',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(*,*) 'case5:',tmp_n
891 ftmp = -hpai*sin(dnei)/dr2/dist
912 ndiff = tmp_n-fnei(i,imin)
915 if (dtmp.ge.15.0d0) then
918 c ntmp = dfaexp( idint(dtmp*1000) + 1 )
922 enei=enei+sccnei(i,imin)*ntmp
924 & sccnei(i,imin)*ntmp*ndiff*2.0d0
925 scc=scc+sccnei(i,imin)
927 c write(*,'(a,1x,2i8,f12.7,i8,3f12.7)')'NEI:',i,imin,tmp_n,
928 c & fnei(i,imin),sccnei(i,imin),enei,scc
931 enei=-enei/scc*snorm_nei*nei_inc*wwnei
932 tmp_fnei=tmp_fnei/scc*snorm_nei*nei_inc*wwnei
934 c if (abs(enei).lt.1.0d-20)then
937 c if (abs(tmp_fnei).lt.1.0d-20) then
946 do j=ishiftca+1,ilastca
947 t2dx(j)=t2dx(j)*tmp_fnei
948 t2dy(j)=t2dy(j)*tmp_fnei
949 t2dz(j)=t2dz(j)*tmp_fnei
952 gdfan(1,n1atom)=gdfan(1,n1atom)+t1dx
953 gdfan(2,n1atom)=gdfan(2,n1atom)+t1dy
954 gdfan(3,n1atom)=gdfan(3,n1atom)+t1dz
956 do j=ishiftca+1,ilastca
957 gdfan(1,j)=gdfan(1,j)+t2dx(j)
958 gdfan(2,j)=gdfan(2,j)+t2dy(j)
959 gdfan(3,j)=gdfan(3,j)+t2dz(j)
972 subroutine edfab(edfabeta)
974 implicit real*8 (a-h,o-z)
977 include 'COMMON.CHAIN'
978 include 'COMMON.DERIV'
982 parameter(PAI=3.14159265358979323846D0)
984 real*8 bx(maxres),by(maxres),bz(maxres)
985 real*8 vbet(maxres,maxres)
986 real*8 shetfx(maxres),shetfy(maxres),shetfz(maxres)
987 real*8 shefx(maxres,12),shefy(maxres,12),shefz(maxres,12)
988 real*8 vbeta,vbetp,vbetm
989 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
991 real*8 dp45,dm45,w_beta
993 common /sheca/ bx,by,bz
994 common /shee/ vbeta,vbet,vbetp,vbetm
995 common /shetf/ shetfx,shetfy,shetfz
996 common /shef/ shefx, shefy, shefz
997 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
999 common /sheconst/ dp45,dm45,w_beta
1000 C End of sheet variables
1003 double precision enebet
1006 bx=0.0d0;by=0.0d0;bz=0.0d0
1007 shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
1011 do i=ishiftca+1,ilastca
1012 bx(i-ishiftca)=c(1,i)
1013 by(i-ishiftca)=c(2,i)
1014 bz(i-ishiftca)=c(3,i)
1021 ULNEX=COS(60.0D0/180.0D0*PAI)
1028 C00=COS((1.0D0+10.0D0/180.0D0)*PAI)
1029 S00=SIN((1.0D0+10.0D0/180.0D0)*PAI)
1035 C END OF INITIALIZATION
1037 nca=ilastca-ishiftca
1039 call angvectors(nca)
1040 call sheetforce(nca,wshet,dfaexp)
1042 c end of sheet energy and force
1045 shetfx(j)=shetfx(j)*beta_inc
1046 shetfy(j)=shetfy(j)*beta_inc
1047 shetfz(j)=shetfz(j)*beta_inc
1048 c write(*,*)'SHETF:',shetfx(j),shetfy(j),shetfz(j)
1051 vbeta=vbeta*beta_inc
1056 gdfab(1,j+ishiftca)=gdfab(1,j+ishiftca)-shetfx(j)
1057 gdfab(2,j+ishiftca)=gdfab(2,j+ishiftca)-shetfy(j)
1058 gdfab(3,j+ishiftca)=gdfab(3,j+ishiftca)-shetfz(j)
1063 C-------------------------------------------------------------------------------
1064 subroutine angvectors(nca)
1065 c implicit real*4(a-h,o-z)
1069 parameter(maxca=800)
1071 parameter(PAI=3.14159265358979323846D0,zero=0.0d0)
1073 real*8 bx(maxca),by(maxca),bz(maxca)
1074 real*8 dis(maxca,maxca)
1075 real*8 apx(maxca),apy(maxca),apz(maxca)
1076 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
1077 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
1078 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
1079 real*8 atx(maxca),aty(maxca),atz(maxca)
1080 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
1081 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
1082 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
1083 real*8 astx(maxca),asty(maxca),astz(maxca)
1084 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1085 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1086 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1088 real*8 cph(maxca),cth(maxca)
1091 integer i, ip, ipp, ip3, j
1092 real*8 rx(maxca, maxca), ry(maxca, maxca), rz(maxca, maxca)
1093 real*8 rix, riy, riz, ripx, ripy, ripz, rippx, rippy, rippz
1094 real*8 gix, giy, giz, gipx, gipy, gipz, gippx, gippy, gippz
1095 real*8 cix, ciy, ciz, cipx, cipy, cipz
1096 real*8 gpcrp_x, gpcrp_y, gpcrp_z, d_gpcrp, gpcrp__g
1097 real*8 d10, d11, d12, d13, d20, d21, d22, d23, d24
1098 real*8 d30, d31, d32, d33, d34, d35, d40, d41, d42, d43
1099 real*8 d_gcr, d_gcr3, d_gmcrim,d_gmcrim3,dgmmcrimm,d_gmmcrimm3
1100 real*8 dg, dg3, dg30, dgm, dgm3, dgmm, dgmm3, dgp, dri
1101 real*8 dri3, drim, drim3, drimm, drip, dripp, g3gmm, g3rim
1102 real*8 g3x, g3y, g3z, d_gmmcrimm, g3rim_,gcr__gm
1103 real*8 gcr_x,gcr_y,gcr_z,ggm,ggp,gmcrim__gmm
1104 real*8 gmcrim_x,gmcrim_y,gmcrim_z,gmmcrimm__gmmm
1105 real*8 gmmcrimm_x,gmmcrimm_y,gmmcrimm_z,gmmgm,gmmr
1106 real*8 gmmx,gmmy,gmmz,gmrp,gmx,gmy,gmz,gpx,gpy,gpz
1107 real*8 grpp,gx,gy,gz
1108 real*8 rim3x,rim3y,rim3z,rimmx,rimmy,rimmz,rimx,rimy,rimz
1109 real*8 sd10,sd11,sd20,sd21,sd22,sd30,sd31,sd32,sd40,sd41
1110 integer inb,nmax,iselect
1112 common /sheca/ bx,by,bz
1113 common /difvec/ rx, ry, rz
1114 common /ulang/ ulcos
1115 common /phys1/ inb,nmax,iselect
1118 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
1119 & apmmz,apm3x,apm3y,apm3z
1120 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
1121 & atmmz,atm3x,atm3y,atm3z
1122 common /coscos/ cph,cth
1123 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1124 & astmmz,astm3x,astm3y,astm3z
1126 C-------------------------------------------------------------------------------
1127 c write(*,*) 'inside angvectors'
1132 cph=zero; cth=zero; sth=zero
1133 apx=zero;apy=zero;apz=zero;apmx=zero;apmy=zero;apmz=zero
1134 apmmx=zero;apmmy=zero;apmmz=zero;apm3x=zero;apm3y=zero;apm3z=zero
1135 atx=zero;aty=zero;atz=zero;atmx=zero;atmy=zero;atmz=zero
1136 atmmx=zero;atmmy=zero;atmmz=zero;atm3x=zero;atm3y=zero;atm3z=zero
1137 astx=zero;asty=zero;astz=zero;astmx=zero;astmy=zero;astmz=zero
1138 astmmx=zero;astmmy=zero;astmmz=zero;astm3x=zero;astm3y=zero
1141 C r[x,y,z] calc and distance calculation
1142 rx=zero;ry=zero;rz=zero
1149 dis(i,j)=sqrt(rx(i,j)**2+ry(i,j)**2+rz(i,j)**2)
1150 c write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
1151 c write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
1152 c write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
1153 c write(*,*) 'dis(i,j):',i,j,dis(i,j)
1156 c end of r[x,y,z] calc
1162 if(dis(i,ip).ge.1.0e-8.and.dis(ip,ipp).ge.1.0e-8) then
1163 ulcos(i)=rx(i,ip)*rx(ip,ipp)+ry(i,ip)*ry(ip,ipp)
1164 $ +rz(i,ip)*rz(ip,ipp)
1165 ulcos(i)=ulcos(i)/(dis(i,ip)*dis(ip,ipp))
1168 c end of virtual bond angle
1169 c write(*,*) 'inside angvectors1'
1180 rippx=bx(ip3)-bx(ipp)
1181 rippy=by(ip3)-by(ipp)
1182 rippz=bz(ip3)-bz(ipp)
1184 gx=riy*ripz-riz*ripy
1185 gy=riz*ripx-rix*ripz
1186 gz=rix*ripy-riy*ripx
1187 gpx=ripy*rippz-ripz*rippy
1188 gpy=ripz*rippx-ripx*rippz
1189 gpz=ripx*rippy-ripy*rippx
1190 gpcrp_x=gpy*ripz-gpz*ripy
1191 gpcrp_y=gpz*ripx-gpx*ripz
1192 gpcrp_z=gpx*ripy-gpy*ripx
1193 d_gpcrp=sqrt(gpcrp_x**2+gpcrp_y**2+gpcrp_z**2)
1194 gpcrp__g=gx*gpy*ripz+gpx*ripy*gz+ripx*gpz*gy
1195 & -gz*gpy*ripx-gpz*ripy*gx-ripz*gpx*gy
1201 gmx=rimy*riz-rimz*riy
1202 gmy=rimz*rix-rimx*riz
1203 gmz=rimx*riy-rimy*rix
1204 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1206 ggm=gmx*gx+gmy*gy+gmz*gz
1207 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1213 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1215 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1216 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1218 c write(*,*) 'inside angvectors2'
1220 rimmx=bx(i-1)-bx(i-2)
1221 rimmy=by(i-1)-by(i-2)
1222 rimmz=bz(i-1)-bz(i-2)
1224 gmmx=rimmy*rimz-rimmz*rimy
1225 gmmy=rimmz*rimx-rimmx*rimz
1226 gmmz=rimmx*rimy-rimmy*rimx
1227 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1229 gmmgm=gmmx*gmx+gmmy*gmy+gmmz*gmz
1230 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1231 gmcrim_x=gmy*rimz-gmz*rimy
1232 gmcrim_y=gmz*rimx-gmx*rimz
1233 gmcrim_z=gmx*rimy-gmy*rimx
1234 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1235 d_gmcrim3=d_gmcrim**3
1236 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1237 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1241 rim3x=bx(i-2)-bx(i-3)
1242 rim3y=by(i-2)-by(i-3)
1243 rim3z=bz(i-2)-bz(i-3)
1244 g3x=rim3y*rimmz-rim3z*rimmy
1245 g3y=rim3z*rimmx-rim3x*rimmz
1246 g3z=rim3x*rimmy-rim3y*rimmx
1247 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1248 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1249 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1250 cc**********************************************************************
1251 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1252 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1253 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1254 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1255 d_gmmcrimm3=d_gmmcrimm**3
1256 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1257 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1264 dg=sqrt(gx**2+gy**2+gz**2)
1265 dgp=sqrt(gpx**2+gpy**2+gpz**2)
1268 ggp=gx*gpx+gy*gpy+gz*gpz
1269 grpp=gx*rippx+gy*rippy+gz*rippz
1271 if(dg.gt.0.0D0.and.dripp.gt.0.0D0.and.dgp.gt.0.0D0
1272 & .and.d_gpcrp.gt.0.0D0) then
1273 cph(i)=grpp/dg/dripp
1275 sth(i)=gpcrp__g/d_gpcrp/dg
1283 c write(*,*) 'inside angvectors3'
1285 if(dgp.gt.0.0D0.and.dg3.gt.0.0D0
1286 & .and.dripp.gt.0.0D0.and.d_gpcrp.gt.0.0D0) then
1289 d12=1.0D0/(dg*dripp)
1290 d13=grpp/(dg3*dripp)
1291 sd10=1.0D0/(d_gpcrp*dg)
1292 sd11=gpcrp__g/(d_gpcrp*dg3)
1302 atx(i)=(ripz*gpy-ripy*gpz)*d10
1303 & -(gy*ripz-gz*ripy)*d11
1304 aty(i)=(ripx*gpz-ripz*gpx)*d10
1305 & -(gz*ripx-gx*ripz)*d11
1306 atz(i)=(ripy*gpx-ripx*gpy)*d10
1307 & -(gx*ripy-gy*ripx)*d11
1308 astx(i)=sd10*(-gpx*ripy**2+ripx*gpz*ripz
1309 & +ripy*gpy*ripx-gpx*ripz**2)
1310 & -sd11*(gy*ripz-gz*ripy)
1311 asty(i)=sd10*(-gpy*ripz**2+gpx*ripy*ripx
1312 & -gpy*ripx**2+gpz*ripy*ripz)
1313 & -sd11*(-gx*ripz+gz*ripx)
1314 astz(i)=sd10*(ripy*gpy*ripz-gpz*ripx**2
1315 & -gpz*ripy**2+ripz*gpx*ripx)
1316 & -sd11*(gx*ripy-gy*ripx)
1317 apx(i)=(ripz*rippy-ripy*rippz)*d12
1318 & -(gy*ripz-gz*ripy)*d13
1319 apy(i)=(ripx*rippz-ripz*rippx)*d12
1320 & -(gz*ripx-gx*ripz)*d13
1321 apz(i)=(ripy*rippx-ripx*rippy)*d12
1322 & -(gx*ripy-gy*ripx)*d13
1334 if(dgm3.gt.0.0D0.and.dg3.gt.0.0D0.and.drip.gt.0.0D0
1335 & .and.d_gcr3.gt.0.0D0) then
1339 d23=1.0D0/(dgm*drip)
1340 d24=gmrp/(dgm3*drip)
1341 sd20=1.0D0/(d_gcr*dgm)
1342 sd21=gcr__gm/(d_gcr3*dgm)
1343 sd22=gcr__gm/(d_gcr*dgm3)
1354 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1355 & -(ciy*gmz-ciz*gmy)*d21
1356 & +(ripy*gz-ripz*gy)*d22
1357 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1358 & -(ciz*gmx-cix*gmz)*d21
1359 & +(ripz*gx-ripx*gz)*d22
1360 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1361 & -(cix*gmy-ciy*gmx)*d21
1362 & +(ripx*gy-ripy*gx)*d22
1363 cc**********************************************************************
1364 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1365 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1366 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1367 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1368 & +gcr_z*(-ripz*rix+gy))
1369 & -sd22*(-gmy*ciz+gmz*ciy)
1371 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1372 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1374 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1375 & -gcr_z*(ripz*riy+gx))
1376 & -sd22*(gmx*ciz-gmz*cix)
1378 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1379 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1381 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1382 & +gcr_z*(ripy*riy+ripx*rix))
1383 & -sd22*(-gmx*ciy+gmy*cix)
1384 cc**********************************************************************
1385 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1386 & -(ciy*gmz-ciz*gmy)*d24
1387 apmy(i)=(ciz*ripx-ripz*cix)*d23
1388 & -(ciz*gmx-cix*gmz)*d24
1389 apmz(i)=(cix*ripy-ripx*ciy)*d23
1390 & -(cix*gmy-ciy*gmx)*d24
1394 if(dgm3.gt.0.0D0.and.dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1395 & .and.d_gmcrim3.gt.0.0D0) then
1396 d30=1.0D0/(dgm*dgmm)
1397 d31=gmmgm/(dgm3*dgmm)
1398 d32=gmmgm/(dgm*dgmm3)
1399 d33=1.0D0/(dgmm*dri)
1400 d34=gmmr/(dgmm3*dri)
1401 d35=gmmr/(dgmm*dri3)
1402 sd30=1.0D0/(d_gmcrim*dgmm)
1403 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1404 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1417 c write(*,*) 'inside angvectors4'
1419 cc**********************************************************************
1420 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1421 & -(ciy*gmz-ciz*gmy)*d31
1422 & -(gmmy*rimmz-gmmz*rimmy)*d32
1423 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1424 & -(ciz*gmx-cix*gmz)*d31
1425 & -(gmmz*rimmx-gmmx*rimmz)*d32
1426 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1427 & -(cix*gmy-ciy*gmx)*d31
1428 & -(gmmx*rimmy-gmmy*rimmx)*d32
1429 cc**********************************************************************
1430 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1431 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1432 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1433 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1434 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1435 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1436 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1438 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1439 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1440 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1441 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1442 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1443 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1444 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1446 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1447 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1448 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1449 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1450 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1451 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1452 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1453 c**********************************************************************
1454 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1455 & -(gmmy*rimmz-gmmz*rimmy)*d34
1457 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1458 & -(gmmz*rimmx-gmmx*rimmz)*d34
1460 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1461 & -(gmmx*rimmy-gmmy*rimmx)*d34
1466 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1467 & .and.drim3.gt.0.0D0
1468 & .and.d_gmmcrimm3.gt.0.0D0) then
1469 d40=1.0D0/(dg30*dgmm)
1470 d41=g3gmm/(dg30*dgmm3)
1471 d42=1.0D0/(dg30*drim)
1472 d43=g3rim_/(dg30*drim3)
1473 sd40=1.0D0/(dg30*d_gmmcrimm)
1474 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1483 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1484 & -(gmmy*rimmz-gmmz*rimmy)*d41
1485 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1486 & -(gmmz*rimmx-gmmx*rimmz)*d41
1487 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1488 & -(gmmx*rimmy-gmmy*rimmx)*d41
1489 cc**********************************************************************
1490 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1491 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1492 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1493 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1495 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1496 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1497 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1498 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1500 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1501 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1502 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1503 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1504 c**********************************************************************
1505 apm3x(i)=g3x*d42-rimx*d43
1506 apm3y(i)=g3y*d42-rimy*d43
1507 apm3z(i)=g3z*d42-rimz*d43
1510 c*******************************************************************************
1512 c write(*,*) 'inside angvectors5'
1517 rimz=bz(i)-bz(i-1)
\r
1518 rimmx=bx(i-1)-bx(i-2)
1519 rimmy=by(i-1)-by(i-2)
1520 rimmz=bz(i-1)-bz(i-2)
1521 rim3x=bx(i-2)-bx(i-3)
1522 rim3y=by(i-2)-by(i-3)
1523 rim3z=bz(i-2)-bz(i-3)
\r
1524 gmmx=rimmy*rimz-rimmz*rimy
1525 gmmy=rimmz*rimx-rimmx*rimz
1526 gmmz=rimmx*rimy-rimmy*rimx
1527 g3x=rim3y*rimmz-rim3z*rimmy
1528 g3y=rim3z*rimmx-rim3x*rimmz
1529 g3z=rim3x*rimmy-rim3y*rimmx
1531 dg30=sqrt(g3x**2+g3y**2+g3z**2)
1532 g3gmm=g3x*gmmx+g3y*gmmy+g3z*gmmz
1533 dgmm=sqrt(gmmx**2+gmmy**2+gmmz**2)
1538 g3rim_=g3x*rimx+g3y*rimy+g3z*rimz
1539 cc**********************************************************************
1540 gmmcrimm_x=gmmy*rimmz-gmmz*rimmy
1541 gmmcrimm_y=gmmz*rimmx-gmmx*rimmz
1542 gmmcrimm_z=gmmx*rimmy-gmmy*rimmx
1543 d_gmmcrimm=sqrt(gmmcrimm_x**2+gmmcrimm_y**2+gmmcrimm_z**2)
1544 d_gmmcrimm3=d_gmmcrimm**3
1545 gmmcrimm__gmmm=g3x*gmmy*rimmz+gmmx*rimmy*g3z+rimmx*gmmz*g3y
1546 & -g3z*gmmy*rimmx-gmmz*rimmy*g3x-rimmz*gmmx*g3y
1548 if(dg30.gt.0.0D0.and.dgmm3.gt.0.0D0
1549 & .and.drim3.gt.0.0D0
1550 & .and.d_gmmcrimm3.gt.0.0D0) then
1551 d40=1.0D0/(dg30*dgmm)
1552 d41=g3gmm/(dg30*dgmm3)
1553 d42=1.0D0/(dg30*drim)
1554 d43=g3rim_/(dg30*drim3)
1555 sd40=1.0D0/(dg30*d_gmmcrimm)
1556 sd41=gmmcrimm__gmmm/(d_gmmcrimm3*dg30)
1565 atm3x(i)=(g3y*rimmz-g3z*rimmy)*d40
1566 & -(gmmy*rimmz-gmmz*rimmy)*d41
1567 atm3y(i)=(g3z*rimmx-g3x*rimmz)*d40
1568 & -(gmmz*rimmx-gmmx*rimmz)*d41
1569 atm3z(i)=(g3x*rimmy-g3y*rimmx)*d40
1570 & -(gmmx*rimmy-gmmy*rimmx)*d41
1571 cc**********************************************************************
1572 astm3x(i)=sd40*(g3x*rimmz**2-rimmx*rimmy*g3y
1573 & -g3z*rimmz*rimmx+rimmy**2*g3x)
1574 & -sd41*(gmmcrimm_x*(rimmz**2+rimmy**2)
1575 & -gmmcrimm_y*rimmy*rimmx-gmmcrimm_z*rimmz*rimmx)
1577 astm3y(i)=sd40*(-rimmz*rimmy*g3z+rimmx**2*g3y
1578 & -rimmx*rimmy*g3x+rimmz**2*g3y)
1579 & -sd41*(-gmmcrimm_x*rimmx*rimmy
1580 & +gmmcrimm_y*(rimmx**2+rimmz**2)-gmmcrimm_z*rimmz*rimmx)
1582 astm3z(i)=sd40*(-g3x*rimmx*rimmz+rimmy**2*g3z
1583 & +g3z*rimmx**2-rimmz*rimmy*g3y)
1584 & -sd41*(-gmmcrimm_x*rimmx*rimmz-gmmcrimm_y*rimmy*rimmz
1585 & +gmmcrimm_z*(rimmy**2+rimmx**2))
1586 cc**********************************************************************
1587 apm3x(i)=g3x*d42-rimx*d43
1588 apm3y(i)=g3y*d42-rimy*d43
1589 apm3z(i)=g3z*d42-rimz*d43
1599 gmx=rimy*riz-rimz*riy
1600 gmy=rimz*rix-rimx*riz
1601 gmz=rimx*riy-rimy*rix
1602 dgm=sqrt(gmx**2+gmy**2+gmz**2)
1606 gmmgm=gmmx*gmx+gmmy*gmy+gmmz+gmz
1607 gmmr=gmmx*rix+gmmy*riy+gmmz*riz
1608 gmcrim_x=gmy*rimz-gmz*rimy
1609 gmcrim_y=gmz*rimx-gmx*rimz
1610 gmcrim_z=gmx*rimy-gmy*rimx
1611 d_gmcrim=sqrt(gmcrim_x**2+gmcrim_y**2+gmcrim_z**2)
1612 d_gmcrim3=d_gmcrim**3
1613 gmcrim__gmm=gmmx*gmy*rimz+gmx*rimy*gmmz+rimx*gmz*gmmy
1614 & -gmmz*gmy*rimx-gmz*rimy*gmmx-rimz*gmx*gmmy
1616 if(dgm3.gt.0.0D0.and.
1617 & dgmm3.gt.0.0D0.and.dri3.gt.0.0D0
1618 & .and.d_gmcrim3.gt.0.0D0) then
1619 d30=1.0D0/(dgm*dgmm)
1620 d31=gmmgm/(dgm3*dgmm)
1621 d32=gmmgm/(dgm*dgmm3)
1622 d33=1.0D0/(dgmm*dri)
1623 d34=gmmr/(dgmm3*dri)
1624 d35=gmmr/(dgmm*dri3)
1625 sd30=1.0D0/(d_gmcrim*dgmm)
1626 sd31=gmcrim__gmm/(d_gmcrim3*dgmm)
1627 sd32=gmcrim__gmm/(d_gmcrim*dgmm3)
1640 cc**********************************************************************
1641 atmmx(i)=(ciy*gmmz-ciz*gmmy-rimmy*gmz+rimmz*gmy)*d30
1642 & -(ciy*gmz-ciz*gmy)*d31
1643 & -(gmmy*rimmz-gmmz*rimmy)*d32
1644 atmmy(i)=(ciz*gmmx-cix*gmmz-rimmz*gmx+rimmx*gmz)*d30
1645 & -(ciz*gmx-cix*gmz)*d31
1646 & -(gmmz*rimmx-gmmx*rimmz)*d32
1647 atmmz(i)=(cix*gmmy-ciy*gmmx-rimmx*gmy+rimmy*gmx)*d30
1648 & -(cix*gmy-ciy*gmx)*d31
1649 & -(gmmx*rimmy-gmmy*rimmx)*d32
1650 cc**********************************************************************
1651 astmmx(i)=sd30*(-gmmx*ciz*rimz-gmx*rimy*rimmy
1652 & +gmz*gmmy+rimx*ciy*gmmy+rimx*gmz*rimmz
1653 & +rimmy*gmy*rimx+gmmz*ciz*rimx-gmmz*gmy
1654 & -ciy*rimy*gmmx-rimz*gmx*rimmz)
1655 & -sd31*(gmcrim_x*(-ciz*rimz-ciy*rimy)
1656 & +gmcrim_y*(ciy*rimx+gmz)+gmcrim_z*(ciz*rimx-gmy))
1657 & -sd32*(gmmy*rimmz-rimmy*gmmz)
1659 astmmy(i)=sd30*(-rimmz*gmy*rimz+ciz*rimy*gmmz
1660 & +gmx*gmmz+gmx*rimy*rimmx-rimx*cix*gmmy
1661 & -rimmx*gmy*rimx+cix*rimy*gmmx-gmz*gmmx
1662 & +gmz*rimy*rimmz-rimz*ciz*gmmy)
1663 & -sd31*(gmcrim_x*(cix*rimy-gmz)
1664 & +gmcrim_y*(-cix*rimx-ciz*rimz)+gmcrim_z*(ciz*rimy+gmx))
1665 & -sd32*(-gmmx*rimmz+rimmx*gmmz)
1667 astmmz(i)=sd30*(rimmy*gmy*rimz+gmmx*cix*rimz
1668 & +gmmx*gmy-ciy*rimy*gmmz-rimx*gmz*rimmx
1669 & -gmmz*cix*rimx-gmz*rimy*rimmy-gmx*gmmy
1670 & +rimz*ciy*gmmy+rimz*gmx*rimmx)
1671 & -sd31*(gmcrim_x*(cix*rimz+gmy)
1672 & +gmcrim_y*(ciy*rimz-gmx)+gmcrim_z*(-ciy*rimy-cix*rimx))
1673 & -sd32*(gmmx*rimmy-rimmx*gmmy)
1674 cc**********************************************************************
1675 apmmx(i)=(riy*rimmz-riz*rimmy-gmmx)*d33
1676 & -(gmmy*rimmz-gmmz*rimmy)*d34
1678 apmmy(i)=(riz*rimmx-rix*rimmz-gmmy)*d33
1679 & -(gmmz*rimmx-gmmx*rimmz)*d34
1681 apmmz(i)=(rix*rimmy-riy*rimmx-gmmz)*d33
1682 & -(gmmx*rimmy-gmmy*rimmx)*d34
1686 c write(*,*) 'inside angvectors6'
1696 gx=riy*ripz-riz*ripy
1697 gy=riz*ripx-rix*ripz
1698 gz=rix*ripy-riy*ripx
1699 ggm=gmx*gx+gmy*gy+gmz*gz
1700 gmrp=gmx*ripx+gmy*ripy+gmz*ripz
1701 dg=sqrt(gx**2+gy**2+gz**2)
1707 d_gcr=sqrt(gcr_x**2+gcr_y**2+gcr_z**2)
1709 gcr__gm=gmx*gy*riz+gx*riy*gmz+rix*gz*gmy
1710 & -gmz*gy*rix-gz*riy*gmx-riz*gx*gmy
1711 if(dgm3.gt.0.0D0.and.
1712 & dg3.gt.0.0D0.and.drip.gt.0.0D0.and.d_gcr3.gt.0.0D0
1717 d23=1.0D0/(dgm*drip)
1718 d24=gmrp/(dgm3*drip)
1719 sd20=1.0D0/(d_gcr*dgm)
1720 sd21=gcr__gm/(d_gcr3*dgm)
1721 sd22=gcr__gm/(d_gcr*dgm3)
1732 atmx(i)=(ciy*gz-ciz*gy-ripy*gmz+ripz*gmy)*d20
1733 & -(ciy*gmz-ciz*gmy)*d21
1734 & +(ripy*gz-ripz*gy)*d22
1735 atmy(i)=(ciz*gx-cix*gz-ripz*gmx+ripx*gmz)*d20
1736 & -(ciz*gmx-cix*gmz)*d21
1737 & +(ripz*gx-ripx*gz)*d22
1738 atmz(i)=(cix*gy-ciy*gx-ripx*gmy+ripy*gmx)*d20
1739 & -(cix*gmy-ciy*gmx)*d21
1740 & +(ripx*gy-ripy*gx)*d22
1741 cc**********************************************************************
1742 astmx(i)=sd20*(gmx*ripz*riz+gx*riy*ciy-gz*gmy
1743 & -rix*ripy*gmy-rix*gz*ciz-ciy*gy*rix-gmz*ripz*rix
1744 & +gmz*gy+ripy*riy*gmx+riz*gx*ciz)
1745 & -sd21*(gcr_x*(ripz*riz+ripy*riy)+gcr_y*(-ripy*rix-gz)
1746 & +gcr_z*(-ripz*rix+gy))
1747 & -sd22*(-gmy*ciz+gmz*ciy)
1749 astmy(i)=sd20*(ciz*gy*riz-ripz*riy*gmz-gx*gmz-gx*riy*cix
1750 & +rix*ripx*gmy+cix*gy*rix-ripx*riy*gmx+gz*gmx-gz*riy*ciz
1752 & -sd21*(gcr_x*(-ripx*riy+gz)+gcr_y*(ripx*rix+ripz*riz)
1753 & -gcr_z*(ripz*riy+gx))
1754 & -sd22*(gmx*ciz-gmz*cix)
1756 astmz(i)=sd20*(-ciy*gy*riz-gmx*ripx*riz-gmx*gy+ripy*riy*gmz
1757 & +rix*gz*cix+gmz*ripx*rix+gz*riy*ciy+gx*gmy-riz*ripy*gmy
1759 & -sd21*(gcr_x*(-ripx*riz-gy)+gcr_y*(-ripy*riz+gx)
1760 & +gcr_z*(ripy*riy+ripx*rix))
1761 & -sd22*(-gmx*ciy+gmy*cix)
1762 cc**********************************************************************
1764 apmx(i)=(ciy*ripz-ripy*ciz)*d23
1765 & -(ciy*gmz-ciz*gmy)*d24
1766 apmy(i)=(ciz*ripx-ripz*cix)*d23
1767 & -(ciz*gmx-cix*gmz)*d24
1768 apmz(i)=(cix*ripy-ripx*ciy)*d23
1769 & -(cix*gmy-ciy*gmx)*d24
1777 c-------------------------------------------------------------------------------
1778 C---------------------------------------------------------------------------------
1779 subroutine sheetforce(nca,wshet,dfaexp)
1782 c this should be matched with dfa.fcm
1784 parameter(maxca=800)
1785 cc**********************************************************************
1788 integer inb,nmax,iselect
1790 real*8 dfaexp(15001)
1792 real*8 vbeta,vbetp,vbetm
1793 real*8 shefx(maxca,12)
1794 real*8 shefy(maxca,12),shefz(maxca,12)
1795 real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
1796 real*8 vbet(maxca,maxca)
1797 real*8 wshet(maxca,maxca)
1798 real*8 bx(maxca),by(maxca),bz(maxca)
1800 common /sheca/ bx,by,bz
1801 common /phys1/ inb,nmax,iselect
1802 common /shef/ shefx,shefy,shefz
1803 common /shee/ vbeta,vbet,vbetp,vbetm
1804 common /shetf/ shetfx,shetfy,shetfz
1821 call sheetene(nca,wshet,dfaexp)
1824 887 format(a,1x,i6,3x,f12.8)
1825 888 format(a,1x,i4,1x,i4,3x,f12.8)
1826 889 format(a,1x,i4,3x,f12.8)
1827 !write(2,*) 'coord : '
1829 !write(2,887) 'bx:',i,bx(i)
1830 !write(2,887) 'by:',i,by(i)
1831 !write(2,887) 'bz:',i,bz(i)
1833 !write(2,*) 'After sheetforce1'
1836 !write(2,888) 'shefx :',i,k,shefx(i,k)
1837 !write(2,888) 'shefy :',i,k,shefy(i,k)
1838 !write(2,888) 'shefz :',i,k,shefz(i,k)
1844 !write(2,*) 'After sheetforce5'
1847 !write(2,888) 'shefx :',i,k,shefx(i,k)
1848 !write(2,888) 'shefy :',i,k,shefy(i,k)
1849 !write(2,888) 'shefz :',i,k,shefz(i,k)
1855 !write(2,*) 'After sheetforce6'
1858 !write(2,888) 'shefx :',i,k,shefx(i,k)
1859 !write(2,888) 'shefy :',i,k,shefy(i,k)
1860 !write(2,888) 'shefz :',i,k,shefz(i,k)
1866 !write(2,*) 'After sheetforce11'
1869 !write(2,888) 'shefx :',i,k,shefx(i,k)
1870 !write(2,888) 'shefy :',i,k,shefy(i,k)
1871 !write(2,888) 'shefz :',i,k,shefz(i,k)
1877 !write(2,*) 'After sheetforce12'
1880 !write(2,888) 'shefx :',i,k,shefx(i,k)
1881 !write(2,888) 'shefy :',i,k,shefy(i,k)
1882 !write(2,888) 'shefz :',i,k,shefz(i,k)
1888 shetfx(i)=shetfx(i)+shefx(i,k)
1889 shetfy(i)=shetfy(i)+shefy(i,k)
1890 shetfz(i)=shetfz(i)+shefz(i,k)
1893 !write(2,*) 'Beta Finished'
1895 !write(2,889) 'shetfx : ',i,shetfx(i)
1896 !write(2,889) 'shetfy : ',i,shetfy(i)
1897 !write(2,889) 'shetfz : ',i,shetfz(i)
1903 c-------------------------------------------------------------------------------
1904 subroutine sheetene(nca,wshet,dfaexp)
1907 parameter(maxca=800)
1908 cc******************************************************************************
1910 real*8 dfaexp(15001)
1911 real*8 dtmp1, dtmp2, dtmp3
1913 real*8 vbet(maxca,maxca)
1914 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
1915 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
1916 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
1917 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
1918 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
1919 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
1920 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
1921 real*8 cph(maxca),cth(maxca)
1922 real*8 rx(maxca,maxca)
1923 real*8 ry(maxca,maxca),rz(maxca,maxca)
1924 real*8 bx(maxca),by(maxca),bz(maxca)
1925 real*8 dis(maxca,maxca)
1927 cc**********************************************************************
1928 real*8 astx(maxca),asty(maxca),astz(maxca)
1929 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
1930 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
1931 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
1933 real*8 wshet(maxca,maxca)
1934 real*8 dp45, dm45, w_beta
1935 real*8 c00, s00, ulnex, dnex, dca,dlhb,ulhb,dshe,dldhb,uldhb
1937 integer i,ip,ipp,j,jp,jpp,inb,nmax,iselect
1939 real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
1941 common /sheca/ bx,by,bz
1942 common /phys1/ inb,nmax,iselect
1944 common /difvec/ rx,ry,rz
1945 common /coscos/ cph,cth
1946 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
1947 & c00,s00,ulnex,dnex
1948 common /sheconst/ dp45,dm45,w_beta
1949 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
1950 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
1951 common /shee/ vbeta,vbet,vbetp,vbetm
1952 common /ulang/ ulcos
1953 cc**********************************************************************
1954 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
1955 & astmmz,astm3x,astm3y,astm3z
1958 real*8 r_pair_mat(maxca,maxca)
1959 common /beta_p/ r_pair_mat
1960 C-------------------------------------------------------------------------------
1964 r_pair_mat(i,j)=wshet(i,j)
1965 c write(*,*) 'r_pair_mat :',i,j,r_pair_mat(i,j)
1980 cc**********************************************************************
1981 y1=(cth(i)*c00+sth(i)*s00-1.0D0)**2
1982 & +(cth(j)*c00+sth(j)*s00-1.0D0)**2
1984 y2=(ulcos(i)-ulnex)**2+(ulcos(ip)-ulnex)**2
1985 & +(ulcos(j)-ulnex)**2+(ulcos(jp)-ulnex)**2
1989 yy1=-0.5d0*(dis(ip,jp)-ulhb)**2/dlhb
1990 yy2=-0.5d0*(dis(ipp,jpp)-ulhb)**2/dlhb
1992 pin1(i,j)=(rx(ip,jp)*rx(ip,ipp)+ry(ip,jp)*ry(ip,ipp)
1993 $ +rz(ip,jp)*rz(ip,ipp))/(dis(ip,jp)*dis(ip,ipp))
1994 pin2(i,j)=(rx(ip,jp)*rx(jp,jpp)+ry(ip,jp)*ry(jp,jpp)
1995 $ +rz(ip,jp)*rz(jp,jpp))/(dis(ip,jp)*dis(jp,jpp))
1996 pin3(i,j)=(rx(ipp,jpp)*rx(ip,ipp)+ry(ipp,jpp)*ry(ip,ipp)
1997 $ +rz(ipp,jpp)*rz(ip,ipp))/(dis(ipp,jpp)*dis(ip,ipp))
1998 pin4(i,j)=(rx(ipp,jpp)*rx(jp,jpp)+ry(ipp,jpp)*ry(jp,jpp)
1999 $ +rz(ipp,jpp)*rz(jp,jpp))/(dis(ipp,jpp)*dis(jp,jpp))
2001 yshe1=pin1(i,j)**2+pin2(i,j)**2
2002 yshe1=-0.5d0*yshe1/dshe
2003 yshe2=pin3(i,j)**2+pin4(i,j)**2
2004 yshe2=-0.5d0*yshe2/dshe
2006 C write(*,*) 'rx(i,j):',i,j,rx(i,j),bx(j),bx(i)
2007 C write(*,*) 'ry(i,j):',i,j,ry(i,j),by(j),by(i)
2008 C write(*,*) 'rz(i,j):',i,j,rz(i,j),bz(j),bz(i)
2009 C write(*,*) 'dis(i,j):',i,j,dis(i,j)
2010 C write(*,*) 'rx(ip,jp):',ip,jp,bx(ip),bx(jp),rx(ip,jp)
2011 C write(*,*) 'rx(ip,ipp):',ip,ipp,bx(ip),bx(ipp),rx(ip,ipp)
2012 C write(*,*) 'pin1:',pin1(i,j)
2013 C write(*,*) 'pin2:',pin2(i,j)
2014 C write(*,*) 'pin3:',pin3(i,j)
2015 C write(*,*) 'pin4:',pin4(i,j)
2018 C write(*,*) 'yy1:',yy1
2019 C write(*,*) 'yy2:',yy2
2020 C write(*,*) 'yshe1:',yshe1
2021 C write(*,*) 'yshe2:',yshe2
2026 dtmp3 = y+yy1+yy2+yshe1+yshe2
2028 C write(*,*)'1', i,j,dtmp1,dtmp2,dtmp3
2029 C write(*,*)'2', y,yy1,yy2
2030 C write(*,*)'3', yshe1,yshe2
2032 if (dtmp3.le.-15.0d0) then
2033 c vbetap(i,j)=-dp45*exp(dtmp3)
2036 c vbetap(i,j)=-dp45*dfaexp(idint(-dtmp3*1000)+1)
2037 vbetap(i,j)=-dp45*exp(dtmp3)
2040 if (dtmp1.le.-15.0d0) then
2041 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2044 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)
2045 c $ *dfaexp(idint(-dtmp1*1000)+1)
2046 vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(dtmp1)
2049 if (dtmp2.le.-15.0d0) then
2050 C vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2053 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)
2054 c $ *dfaexp(idint(-dtmp2*1000)+1)
2055 vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(dtmp2)
2058 c vbetap(i,j)=-dp45*exp(y+yy1+yy2+yshe1+yshe2)
2059 c vbetap1(i,j)=-r_pair_mat(i+1,j+1)*exp(y+yy1+yshe1)
2060 c vbetap2(i,j)=-r_pair_mat(i+2,j+2)*exp(y+yy2+yshe2)
2062 ! write(*,*) 'r_pair_mat>',i+1,j+1,r_pair_mat(i+1,j+1)
2063 ! write(*,*) 'r_pair_mat>',i+2,j+2,r_pair_mat(i+2,j+2)
2066 yy1=-0.5d0*(dis(ip,jpp)-ulhb)**2/dlhb
2067 yy2=-0.5d0*(dis(ipp,jp)-ulhb)**2/dlhb
2069 pina1(i,j)=(rx(ip,jpp)*rx(ip,ipp)+ry(ip,jpp)*ry(ip,ipp)
2070 $ +rz(ip,jpp)*rz(ip,ipp))/(dis(ip,jpp)*dis(ip,ipp))
2071 pina2(i,j)=(rx(ip,jpp)*rx(jp,jpp)+ry(ip,jpp)*ry(jp,jpp)
2072 $ +rz(ip,jpp)*rz(jp,jpp))/(dis(ip,jpp)*dis(jp,jpp))
2073 pina3(i,j)=(rx(jp,ipp)*rx(ip,ipp)+ry(jp,ipp)*ry(ip,ipp)
2074 $ +rz(jp,ipp)*rz(ip,ipp))/(dis(jp,ipp)*dis(ip,ipp))
2075 pina4(i,j)=(rx(jp,ipp)*rx(jp,jpp)+ry(jp,ipp)*ry(jp,jpp)
2076 $ +rz(jp,ipp)*rz(jp,jpp))/(dis(jp,ipp)*dis(jp,jpp))
2078 yshe1=pina1(i,j)**2+pina2(i,j)**2
2079 yshe1=-0.5d0*yshe1/dshe
2080 yshe2=pina3(i,j)**2+pina4(i,j)**2
2081 yshe2=-0.5d0*yshe2/dshe
2083 C write(*,*) 'pina1:',pina1(i,j)
2084 C write(*,*) 'pina2:',pina2(i,j)
2085 C write(*,*) 'pina3:',pina3(i,j)
2086 C write(*,*) 'pina4:',pina4(i,j)
2087 C write(*,*) 'yshe1:',yshe1
2088 C write(*,*) 'yshe2:',yshe2
2089 C write(*,*) 'dshe:',dshe
2091 dtmp3=y+yy1+yy2+yshe1+yshe2
2095 if(dtmp3 .le. -15.0d0) then
2096 c vbetam(i,j)=-dm45*exp(dtmp3)
2099 c vbetam(i,j)=-dm45*dfaexp(idint(-dtmp3*1000)+1)
2100 vbetam(i,j)=-dm45*exp(dtmp3)
2103 if(dtmp1 .le. -15.0d0) then
2104 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2107 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)
2108 c $ *dfaexp(idint(-dtmp1*1000)+1)
2109 vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(dtmp1)
2112 if(dtmp2.le.-15.0d0) then
2113 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2116 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)
2117 c $ *dfaexp(idint(-dtmp2*1000)+1)
2118 vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(dtmp2)
2121 c vbetam(i,j)=-dm45*exp(y+yy1+yy2+yshe1+yshe2)
2122 c vbetam1(i,j)=-r_pair_mat(i+1,j+2)*exp(y+yy1+yshe1)
2123 c vbetam2(i,j)=-r_pair_mat(i+2,j+1)*exp(y+yy2+yshe2)
2125 ! write(*,*) 'r_pair_mat>',i+1,j+2,r_pair_mat(i+1,j+2)
2126 ! write(*,*) 'r_pair_mat>',i+2,j+1,r_pair_mat(i+2,j+1)
2128 uup = vbetap(i,j)+vbetap1(i,j)+vbetap2(i,j)
2129 uum = vbetam(i,j)+vbetam1(i,j)+vbetam2(i,j)
2131 c write(*,*) 'uup,uum:', uup, uum
2133 c uup=vbetap1(i,j)+vbetap2(i,j)
2134 c uum=vbetam1(i,j)+vbetam2(i,j)
2139 vbeta=vbeta+vbet(i,j)
2141 c write(*,*) 'uup,uum:',uup,uum
2142 c write(*,*) 'vbetap(i,j):',vbetap(i,j)
2143 c write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2144 c write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2145 c write(*,*) 'vbetam(i,j):',vbetam(i,j)
2146 c write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2147 c write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2148 c write(*,*) 'uup:',uup
2149 c write(*,*) 'uum:',uum
2150 c write(*,*) 'vbetp:',vbetp
2151 c write(*,*) 'vbetm:',vbetm
2152 c write(*,*) 'vbet(i,j):',vbet(i,j)
2160 ! write(*,*) 'I,J:', i,j
2161 ! write(*,*) 'vbetap(i,j):',vbetap(i,j)
2162 ! write(*,*) 'vbetap1(i,j):',vbetap1(i,j)
2163 ! write(*,*) 'vbetap2(i,j):',vbetap2(i,j)
2164 ! write(*,*) 'vbetam(i,j):',vbetam(i,j)
2165 ! write(*,*) 'vbetam1(i,j):',vbetam1(i,j)
2166 ! write(*,*) 'vbetam2(i,j):',vbetam2(i,j)
2167 ! write(*,*) 'vbet(i,j):',vbet(i,j)
2173 c-------------------------------------------------------------------------------
2174 subroutine sheetforce1
2177 parameter(maxca=800)
2178 cc**********************************************************************
2179 real*8 vbet(maxca,maxca)
2180 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2181 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2182 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2183 real*8 cph(maxca),cth(maxca)
2184 real*8 rx(maxca,maxca)
2185 real*8 ry(maxca,maxca),rz(maxca,maxca)
2186 real*8 bx(maxca),by(maxca),bz(maxca)
2187 real*8 dis(maxca,maxca)
2188 real*8 shefx(maxca,12)
2189 real*8 shefy(maxca,12),shefz(maxca,12)
2190 real*8 atx(maxca),aty(maxca),atz(maxca)
2191 real*8 atmx(maxca),atmy(maxca),atmz(maxca)
2192 real*8 atmmx(maxca),atmmy(maxca),atmmz(maxca)
2193 real*8 atm3x(maxca),atm3y(maxca),atm3z(maxca)
2194 real*8 apx(maxca),apy(maxca),apz(maxca)
2195 real*8 apmx(maxca),apmy(maxca),apmz(maxca)
2196 real*8 apmmx(maxca),apmmy(maxca),apmmz(maxca)
2197 real*8 apm3x(maxca),apm3y(maxca),apm3z(maxca)
2199 real*8 astx(maxca),asty(maxca),astz(maxca)
2200 real*8 astmx(maxca),astmy(maxca),astmz(maxca)
2201 real*8 astmmx(maxca),astmmy(maxca),astmmz(maxca)
2202 real*8 astm3x(maxca),astm3y(maxca),astm3z(maxca)
2204 real*8 w_beta,dp45, dm45
2205 real*8 vbeta, vbetp, vbetm
2206 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2207 $ c00,s00,ulnex,dnex
2208 integer inb,nmax,iselect
2210 common /phys1/ inb,nmax,iselect
2212 common /difvec/ rx,ry,rz
2213 common /coscos/ cph,cth
2214 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2215 $ c00,s00,ulnex,dnex
2216 common /sheconst/ dp45,dm45,w_beta
2217 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2218 common /angvt/ atx,aty,atz,atmx,atmy,atmz,atmmx,atmmy,
2219 $ atmmz,atm3x,atm3y,atm3z
2220 common /angvp/ apx,apy,apz,apmx,apmy,apmz,apmmx,apmmy,
2221 $ apmmz,apm3x,apm3y,apm3z
2222 common /shef/ shefx,shefy,shefz
2223 common /shee/ vbeta,vbet,vbetp,vbetm
2224 common /ulang/ ulcos
2225 c c**********************************************************************
2226 common /angvt2/ astx,asty,astz,astmx,astmy,astmz,astmmx,astmmy,
2227 $ astmmz,astm3x,astm3y,astm3z
2229 C--------------------------------------------------------------------------------
2231 integer i,j,im3,imm,im,ip,ipp,jm,jmm,jm3,jp,jpp
2232 real*8 c1,v1,cc1,dmm,dmm__,fx,fy,fz,c2,v2,dmm1
2233 real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
2234 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
2235 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
2236 C--------------------------------------------------------------------------------
2241 c1=(cth(im3)*c00+sth(im3)*s00-1)/dca
2246 cc1=(ulcos(imm)-ulnex)/dnex
2247 dmm=cc1/(dis(imm,im)*dis(im,i))
2248 dmm__=cc1*ulcos(imm)/dis(im,i)**2
2249 fx=rx(imm,im)*dmm-rx(im,i)*dmm__
2250 fy=ry(imm,im)*dmm-ry(im,i)*dmm__
2251 fz=rz(imm,im)*dmm-rz(im,i)*dmm__
2252 fx=fx+(atm3x(i)*c00+astm3x(i)*s00)*c1
2253 fy=fy+(atm3y(i)*c00+astm3y(i)*s00)*c1
2254 fz=fz+(atm3z(i)*c00+astm3z(i)*s00)*c1
2264 c2=(cth(imm)*c00+sth(imm)*s00-1)/dca
2269 cc1=(ulcos(imm)-ulnex)/dnex
2270 cc2=(ulcos(im)-ulnex)/dnex
2271 dmm1=cc1/(dis(imm,im)*dis(im,i))
2272 dmm2=cc2/(dis(im,i)*dis(i,ip))
2273 dmm1__=cc1*ulcos(imm)/dis(im,i)**2
2274 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2275 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2276 cc**********************************************************************
2277 fx=rx(imm,im)*dmm1-rx(im,i)*dmm1__+rx(i,ip)*dmm2-rx(im,i)*dmm2
2278 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2
2279 fy=ry(imm,im)*dmm1-ry(im,i)*dmm1__+ry(i,ip)*dmm2-ry(im,i)*dmm2
2280 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2
2281 fz=rz(imm,im)*dmm1-rz(im,i)*dmm1__+rz(i,ip)*dmm2-rz(im,i)*dmm2
2282 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2
2283 fx=fx+(atmmx(i)*c00+astmmx(i)*s00)*c2
2284 fy=fy+(atmmy(i)*c00+astmmy(i)*s00)*c2
2285 fz=fz+(atmmz(i)*c00+astmmz(i)*s00)*c2
2294 c3=(cth(im)*c00+sth(im)*s00-1)/dca
2299 cc2=(ulcos(im)-ulnex)/dnex
2300 cc3=(ulcos(i)-ulnex)/dnex
2301 dmm2=cc2/(dis(im,i)*dis(i,ip))
2302 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2303 dmm2_1=cc2*ulcos(im)/dis(im,i)**2
2304 dmm2_2=cc2*ulcos(im)/dis(i,ip)**2
2305 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2306 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm2-rx(im,i)*dmm2
2307 $ -rx(im,i)*dmm2_1+rx(i,ip)*dmm2_2+rx(i,ip)*dmm3__
2308 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm2-ry(im,i)*dmm2
2309 $ -ry(im,i)*dmm2_1+ry(i,ip)*dmm2_2+ry(i,ip)*dmm3__
2310 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm2-rz(im,i)*dmm2
2311 $ -rz(im,i)*dmm2_1+rz(i,ip)*dmm2_2+rz(i,ip)*dmm3__
2312 fx=fx+(atmx(i)*c00+astmx(i)*s00)*c3
2313 fy=fy+(atmy(i)*c00+astmy(i)*s00)*c3
2314 fz=fz+(atmz(i)*c00+astmz(i)*s00)*c3
2322 c4=(cth(i)*c00+sth(i)*s00-1)/dca
2327 cc3=(ulcos(i)-ulnex)/dnex
2328 dmm3=cc3/(dis(i,ip)*dis(ip,ipp))
2329 dmm3__=cc3*ulcos(i)/dis(i,ip)**2
2330 fx=-rx(ip,ipp)*dmm3+rx(i,ip)*dmm3__
2331 fy=-ry(ip,ipp)*dmm3+ry(i,ip)*dmm3__
2332 fz=-rz(ip,ipp)*dmm3+rz(i,ip)*dmm3__
2333 fx=fx+(atx(i)*c00+astx(i)*s00)*c4
2334 fy=fy+(aty(i)*c00+asty(i)*s00)*c4
2335 fz=fz+(atz(i)*c00+astz(i)*s00)*c4
2344 c7=(cth(jm3)*c00+sth(jm3)*s00-1)/dca
2349 cc7=(ulcos(jmm)-ulnex)/dnex
2350 dmm=cc7/(dis(jmm,jm)*dis(jm,j))
2351 dmm__=cc7*ulcos(jmm)/dis(jm,j)**2
2352 fx=rx(jmm,jm)*dmm-rx(jm,j)*dmm__
2353 fy=ry(jmm,jm)*dmm-ry(jm,j)*dmm__
2354 fz=rz(jmm,jm)*dmm-rz(jm,j)*dmm__
2355 fx=fx+(atm3x(j)*c00+astm3x(j)*s00)*c7
2356 fy=fy+(atm3y(j)*c00+astm3y(j)*s00)*c7
2357 fz=fz+(atm3z(j)*c00+astm3z(j)*s00)*c7
2366 c8=(cth(jmm)*c00+sth(jmm)*s00-1)/dca
2371 cc7=(ulcos(jmm)-ulnex)/dnex
2372 cc8=(ulcos(jm)-ulnex)/dnex
2373 dmm7=cc7/(dis(jmm,jm)*dis(jm,j))
2374 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2375 dmm7__=cc7*ulcos(jmm)/dis(jm,j)**2
2376 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2377 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2378 fx=rx(jmm,jm)*dmm7+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2379 $ -rx(jm,j)*dmm7__-rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2
2380 fy=ry(jmm,jm)*dmm7+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2381 $ -ry(jm,j)*dmm7__-ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2
2382 fz=rz(jmm,jm)*dmm7+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2383 $ -rz(jm,j)*dmm7__-rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2
2384 fx=fx+(atmmx(j)*c00+astmmx(j)*s00)*c8
2385 fy=fy+(atmmy(j)*c00+astmmy(j)*s00)*c8
2386 fz=fz+(atmmz(j)*c00+astmmz(j)*s00)*c8
2396 c9=(cth(jm)*c00+sth(jm)*s00-1)/dca
2401 cc8=(ulcos(jm)-ulnex)/dnex
2402 cc9=(ulcos(j)-ulnex)/dnex
2403 dmm8=cc8/(dis(jm,j)*dis(j,jp))
2404 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2405 dmm8_1=cc8*ulcos(jm)/dis(jm,j)**2
2406 dmm8_2=cc8*ulcos(jm)/dis(j,jp)**2
2407 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2408 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm8-rx(jm,j)*dmm8
2409 $ -rx(jm,j)*dmm8_1+rx(j,jp)*dmm8_2+rx(j,jp)*dmm9__
2410 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm8-ry(jm,j)*dmm8
2411 $ -ry(jm,j)*dmm8_1+ry(j,jp)*dmm8_2+ry(j,jp)*dmm9__
2412 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm8-rz(jm,j)*dmm8
2413 $ -rz(jm,j)*dmm8_1+rz(j,jp)*dmm8_2+rz(j,jp)*dmm9__
2414 fx=fx+(atmx(j)*c00+astmx(j)*s00)*c9
2415 fy=fy+(atmy(j)*c00+astmy(j)*s00)*c9
2416 fz=fz+(atmz(j)*c00+astmz(j)*s00)*c9
2425 c10=(cth(j)*c00+sth(j)*s00-1)/dca
2430 cc9=(ulcos(j)-ulnex)/dnex
2431 dmm9=cc9/(dis(j,jp)*dis(jp,jpp))
2432 dmm9__=cc9*ulcos(j)/dis(j,jp)**2
2433 fx=-rx(jp,jpp)*dmm9+rx(j,jp)*dmm9__
2434 fy=-ry(jp,jpp)*dmm9+ry(j,jp)*dmm9__
2435 fz=-rz(jp,jpp)*dmm9+rz(j,jp)*dmm9__
2436 fx=fx+(atx(j)*c00+astx(j)*s00)*c10
2437 fy=fy+(aty(j)*c00+asty(j)*s00)*c10
2438 fz=fz+(atz(j)*c00+astz(j)*s00)*c10
2446 c----------------------------------------------------------------------------
2447 subroutine sheetforce5
2450 parameter(maxca=800)
2451 cc**********************************************************************
2452 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2453 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2454 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2455 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2456 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2457 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2458 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2459 real*8 rx(maxca,maxca)
2460 real*8 ry(maxca,maxca),rz(maxca,maxca)
2461 real*8 bx(maxca),by(maxca),bz(maxca)
2462 real*8 dis(maxca,maxca)
2463 real*8 shefx(maxca,12),shefy(maxca,12)
2464 real*8 shefz(maxca,12)
2465 real*8 dp45,dm45,w_beta
2466 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2467 $ c00,s00,ulnex,dnex
2468 integer inb,nmax,iselect
2469 cc**********************************************************************
2470 common /phys1/ inb,nmax,iselect
2472 common /difvec/ rx,ry,rz
2473 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2474 $ c00,s00,ulnex,dnex
2475 common /sheconst/ dp45,dm45,w_beta
2476 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2477 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2478 common /shef/ shefx,shefy,shefz
2479 c********************************************************************************
2481 integer i,imm,im,jp,jpp,j
2482 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2483 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2484 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z
2485 real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
2486 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2487 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
2488 c********************************************************************************
2496 yy1=-(dis(i,jpp)-ulhb)/dlhb
2497 y1x=rx(jpp,i)/dis(i,jpp)
2498 y1y=ry(jpp,i)/dis(i,jpp)
2499 y1z=rz(jpp,i)/dis(i,jpp)
2504 yy33=1.0D0/(dis(im,jp)*dis(im,i))
2505 yyy3=pin1(imm,j)/(dis(im,i)**2)
2506 yy3=-pin1(imm,j)/dshe
2507 y3x=(yy33*rx(im,jp)-yyy3*rx(im,i))*yy3
2508 y3y=(yy33*ry(im,jp)-yyy3*ry(im,i))*yy3
2509 y3z=(yy33*rz(im,jp)-yyy3*rz(im,i))*yy3
2511 yy44=1.0D0/(dis(i,jpp)*dis(im,i))
2512 yyy4a=pin3(imm,j)/(dis(i,jpp)**2)
2513 yyy4b=pin3(imm,j)/(dis(im,i)**2)
2514 yy4=-pin3(imm,j)/dshe
2515 y4x=(yy44*(rx(i,jpp)-rx(im,i))+yyy4a*rx(i,jpp)
2516 $ -yyy4b*rx(im,i))*yy4
2517 y4y=(yy44*(ry(i,jpp)-ry(im,i))+yyy4a*ry(i,jpp)
2518 $ -yyy4b*ry(im,i))*yy4
2519 y4z=(yy44*(rz(i,jpp)-rz(im,i))+yyy4a*rz(i,jpp)
2520 $ -yyy4b*rz(im,i))*yy4
2523 yy55=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2524 yyy5=pin4(imm,j)/(dis(i,jpp)**2)
2525 yy5=-pin4(imm,j)/dshe
2526 y5x=(-yy55*rx(jp,jpp)+yyy5*rx(i,jpp))*yy5
2527 y5y=(-yy55*ry(jp,jpp)+yyy5*ry(i,jpp))*yy5
2528 y5z=(-yy55*rz(jp,jpp)+yyy5*rz(i,jpp))*yy5
2541 shefx(i,5)=shefx(i,5)-sx*vbetap(imm,j)
2542 $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2543 shefy(i,5)=shefy(i,5)-sy*vbetap(imm,j)
2544 $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2545 shefz(i,5)=shefz(i,5)-sz*vbetap(imm,j)
2546 $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2548 ! shefx(i,5)=shefx(i,5)
2549 ! $ -sx1*vbetap1(imm,j)-sx2*vbetap2(imm,j)
2550 ! shefy(i,5)=shefy(i,5)
2551 ! $ -sy1*vbetap1(imm,j)-sy2*vbetap2(imm,j)
2552 ! shefz(i,5)=shefz(i,5)
2553 ! $ -sz1*vbetap1(imm,j)-sz2*vbetap2(imm,j)
2555 yy6=-(dis(i,jp)-uldhb)/dldhb
2556 y6x=rx(jp,i)/dis(i,jp)
2557 y6y=ry(jp,i)/dis(i,jp)
2558 y6z=rz(jp,i)/dis(i,jp)
2563 yy88=1.0D0/(dis(im,jpp)*dis(im,i))
2564 yyy8=pina1(imm,j)/(dis(im,i)**2)
2565 yy8=-pina1(imm,j)/dshe
2566 y8x=(yy88*rx(im,jpp)-yyy8*rx(im,i))*yy8
2567 y8y=(yy88*ry(im,jpp)-yyy8*ry(im,i))*yy8
2568 y8z=(yy88*rz(im,jpp)-yyy8*rz(im,i))*yy8
2570 yy99=1.0D0/(dis(jp,i)*dis(im,i))
2571 yyy9a=pina3(imm,j)/(dis(jp,i)**2)
2572 yyy9b=pina3(imm,j)/(dis(im,i)**2)
2573 yy9=-pina3(imm,j)/dshe
2574 y9x=(yy99*(rx(jp,i)+rx(im,i))-yyy9a*rx(jp,i)
2575 $ -yyy9b*rx(im,i))*yy9
2576 y9y=(yy99*(ry(jp,i)+ry(im,i))-yyy9a*ry(jp,i)
2577 $ -yyy9b*ry(im,i))*yy9
2578 y9z=(yy99*(rz(jp,i)+rz(im,i))-yyy9a*rz(jp,i)
2579 $ -yyy9b*rz(im,i))*yy9
2581 yy1010=1.0D0/(dis(jp,i)*dis(jp,jpp))
2582 yyy10=pina4(imm,j)/(dis(jp,i)**2)
2583 yy10=-pina4(imm,j)/dshe
2584 y10x=(yy1010*rx(jp,jpp)-yyy10*rx(jp,i))*yy10
2585 y10y=(yy1010*ry(jp,jpp)-yyy10*ry(jp,i))*yy10
2586 y10z=(yy1010*rz(jp,jpp)-yyy10*rz(jp,i))*yy10
2588 sx=y66x+y8x+y9x+y10x
2589 sy=y66y+y8y+y9y+y10y
2590 sz=y66z+y8z+y9z+y10z
2599 shefx(i,5)=shefx(i,5)-sx*vbetam(imm,j)
2600 $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2601 shefy(i,5)=shefy(i,5)-sy*vbetam(imm,j)
2602 $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2603 shefz(i,5)=shefz(i,5)-sz*vbetam(imm,j)
2604 $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2606 ! shefx(i,5)=shefx(i,5)
2607 ! $ -sx1*vbetam1(imm,j)-sx2*vbetam2(imm,j)
2608 ! shefy(i,5)=shefy(i,5)
2609 ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
2610 ! shefz(i,5)=shefz(i,5)
2611 ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
2618 c--------------------------------------------------------------------------c
2619 subroutine sheetforce6
2622 parameter(maxca=800)
2623 cc**********************************************************************
2624 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2625 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2626 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2627 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2628 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2629 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2630 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2631 real*8 rx(maxca,maxca)
2632 real*8 ry(maxca,maxca),rz(maxca,maxca)
2633 real*8 bx(maxca),by(maxca),bz(maxca)
2634 real*8 dis(maxca,maxca)
2635 real*8 shefx(maxca,12),shefy(maxca,12)
2636 real*8 shefz(maxca,12)
2637 real*8 dp45,dm45,w_beta
2638 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2639 $ c00,s00,ulnex,dnex
2640 integer inb,nmax,iselect
2641 cc**********************************************************************
2642 common /phys1/ inb,nmax,iselect
2644 common /difvec/ rx,ry,rz
2645 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2646 $ c00,s00,ulnex,dnex
2647 common /sheconst/ dp45,dm45,w_beta
2648 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2649 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2650 common /shef/ shefx,shefy,shefz
2651 cc**********************************************************************
2653 integer i,imm,im,jp,jpp,j,ip
2654 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2655 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2656 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
2657 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
2658 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
2659 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
2660 C********************************************************************************
2668 yy1=-(dis(i,jp)-ulhb)/dlhb
2669 y1x=rx(jp,i)/dis(i,jp)
2670 y1y=ry(jp,i)/dis(i,jp)
2671 y1z=rz(jp,i)/dis(i,jp)
2676 yy33=1.0D0/(dis(i,jp)*dis(i,ip))
2677 yyy3a=pin1(im,j)/(dis(i,jp)**2)
2678 yyy3b=pin1(im,j)/(dis(i,ip)**2)
2679 yy3=-pin1(im,j)/dshe
2680 y3x=(-yy33*(rx(i,ip)+rx(i,jp))+yyy3a*rx(i,jp)
2681 $ +yyy3b*rx(i,ip))*yy3
2682 y3y=(-yy33*(ry(i,ip)+ry(i,jp))+yyy3a*ry(i,jp)
2683 $ +yyy3b*ry(i,ip))*yy3
2684 y3z=(-yy33*(rz(i,ip)+rz(i,jp))+yyy3a*rz(i,jp)
2685 $ +yyy3b*rz(i,ip))*yy3
2687 yy44=1.0D0/(dis(i,jp)*dis(jp,jpp))
2688 yyy4=pin2(im,j)/(dis(i,jp)**2)
2689 yy4=-pin2(im,j)/dshe
2690 y4x=(-yy44*rx(jp,jpp)+yyy4*rx(i,jp))*yy4
2691 y4y=(-yy44*ry(jp,jpp)+yyy4*ry(i,jp))*yy4
2692 y4z=(-yy44*rz(jp,jpp)+yyy4*rz(i,jp))*yy4
2694 yy55=1.0D0/(dis(ip,jpp)*dis(i,ip))
2695 yyy5=pin3(im,j)/(dis(i,ip)**2)
2696 yy5=-pin3(im,j)/dshe
2697 y5x=(-yy55*rx(ip,jpp)+yyy5*rx(i,ip))*yy5
2698 y5y=(-yy55*ry(ip,jpp)+yyy5*ry(i,ip))*yy5
2699 y5z=(-yy55*rz(ip,jpp)+yyy5*rz(i,ip))*yy5
2712 shefx(i,6)=shefx(i,6)-sx*vbetap(im,j)
2713 $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
2714 shefy(i,6)=shefy(i,6)-sy*vbetap(im,j)
2715 $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
2716 shefz(i,6)=shefz(i,6)-sz*vbetap(im,j)
2717 $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
2718 ! shefx(i,6)=shefx(i,6)
2719 ! $ -sx1*vbetap1(im,j)-sx2*vbetap2(im,j)
2720 ! shefy(i,6)=shefy(i,6)
2721 ! $ -sy1*vbetap1(im,j)-sy2*vbetap2(im,j)
2722 ! shefz(i,6)=shefz(i,6)
2723 ! $ -sz1*vbetap1(im,j)-sz2*vbetap2(im,j)
2725 yy6=-(dis(jpp,i)-uldhb)/dldhb
2726 y6x=rx(jpp,i)/dis(jpp,i)
2727 y6y=ry(jpp,i)/dis(jpp,i)
2728 y6z=rz(jpp,i)/dis(jpp,i)
2733 yy88=1.0D0/(dis(i,jpp)*dis(i,ip))
2734 yyy8a=pina1(im,j)/(dis(i,jpp)**2)
2735 yyy8b=pina1(im,j)/(dis(i,ip)**2)
2736 yy8=-pina1(im,j)/dshe
2737 y8x=(-yy88*(rx(i,jpp)+rx(i,ip))+yyy8a*rx(i,jpp)
2738 $ +yyy8b*rx(i,ip))*yy8
2739 y8y=(-yy88*(ry(i,jpp)+ry(i,ip))+yyy8a*ry(i,jpp)
2740 $ +yyy8b*ry(i,ip))*yy8
2741 y8z=(-yy88*(rz(i,jpp)+rz(i,ip))+yyy8a*rz(i,jpp)
2742 $ +yyy8b*rz(i,ip))*yy8
2744 yy99=1.0D0/(dis(i,jpp)*dis(jp,jpp))
2745 yyy9=pina2(im,j)/(dis(i,jpp)**2)
2746 yy9=-pina2(im,j)/dshe
2747 y9x=(-yy99*rx(jp,jpp)+yyy9*rx(i,jpp))*yy9
2748 y9y=(-yy99*ry(jp,jpp)+yyy9*ry(i,jpp))*yy9
2749 y9z=(-yy99*rz(jp,jpp)+yyy9*rz(i,jpp))*yy9
2751 yy1010=1.0D0/(dis(jp,ip)*dis(i,ip))
2752 yyy10=pina3(im,j)/(dis(i,ip)**2)
2753 yy10=-pina3(im,j)/dshe
2754 y10x=(-yy1010*rx(jp,ip)+yyy10*rx(i,ip))*yy10
2755 y10y=(-yy1010*ry(jp,ip)+yyy10*ry(i,ip))*yy10
2756 y10z=(-yy1010*rz(jp,ip)+yyy10*rz(i,ip))*yy10
2758 sx=y66x+y8x+y9x+y10x
2759 sy=y66y+y8y+y9y+y10y
2760 sz=y66z+y8z+y9z+y10z
2769 shefx(i,6)=shefx(i,6)-sx*vbetam(im,j)
2770 $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
2771 shefy(i,6)=shefy(i,6)-sy*vbetam(im,j)
2772 $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
2773 shefz(i,6)=shefz(i,6)-sz*vbetam(im,j)
2774 $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
2776 ! shefx(i,6)=shefx(i,6)
2777 ! $ -sx1*vbetam1(im,j)-sx2*vbetam2(im,j)
2778 ! shefy(i,6)=shefy(i,6)
2779 ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
2780 ! shefz(i,6)=shefz(i,6)
2781 ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
2788 c-----------------------------------------------------------------------
2789 subroutine sheetforce11
2792 parameter(maxca=800)
2793 cc**********************************************************************
2794 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2795 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2796 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2797 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2798 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2799 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2800 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2801 real*8 rx(maxca,maxca)
2802 real*8 ry(maxca,maxca),rz(maxca,maxca)
2803 real*8 bx(maxca),by(maxca),bz(maxca)
2804 real*8 dis(maxca,maxca)
2805 real*8 shefx(maxca,12),shefy(maxca,12)
2806 real*8 shefz(maxca,12)
2807 real*8 dp45,dm45,w_beta
2808 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2809 $ c00,s00,ulnex,dnex
2810 integer inb,nmax,iselect
2811 cc**********************************************************************
2812 common /phys1/ inb,nmax,iselect
2814 common /difvec/ rx,ry,rz
2815 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2816 $ c00,s00,ulnex,dnex
2817 common /sheconst/ dp45,dm45,w_beta
2818 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2819 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2820 common /shef/ shefx,shefy,shefz
2821 C********************************************************************************
2823 integer j,jm,jmm,ip,i,ipp
2824 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2825 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y
2826 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y
2827 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
2828 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
2829 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
2830 C********************************************************************************
2839 yy1=-(dis(ipp,j)-ulhb)/dlhb
2840 y1x=rx(ipp,j)/dis(ipp,j)
2841 y1y=ry(ipp,j)/dis(ipp,j)
2842 y1z=rz(ipp,j)/dis(ipp,j)
2847 yy33=1.0D0/(dis(ip,jm)*dis(jm,j))
2848 yyy3=pin2(i,jmm)/(dis(jm,j)**2)
2849 yy3=-pin2(i,jmm)/dshe
2850 y3x=(yy33*rx(ip,jm)-yyy3*rx(jm,j))*yy3
2851 y3y=(yy33*ry(ip,jm)-yyy3*ry(jm,j))*yy3
2852 y3z=(yy33*rz(ip,jm)-yyy3*rz(jm,j))*yy3
2854 yy44=1.0D0/(dis(ipp,j)*dis(ip,ipp))
2855 yyy4=pin3(i,jmm)/(dis(ipp,j)**2)
2856 yy4=-pin3(i,jmm)/dshe
2857 y4x=(yy44*rx(ip,ipp)-yyy4*rx(ipp,j))*yy4
2858 y4y=(yy44*ry(ip,ipp)-yyy4*ry(ipp,j))*yy4
2859 y4z=(yy44*rz(ip,ipp)-yyy4*rz(ipp,j))*yy4
2861 yy55=1.0D0/(dis(ipp,j)*dis(jm,j))
2862 yyy5a=pin4(i,jmm)/(dis(ipp,j)**2)
2863 yyy5b=pin4(i,jmm)/(dis(jm,j)**2)
2864 yy5=-pin4(i,jmm)/dshe
2865 y5x=(yy55*(rx(jm,j)+rx(ipp,j))-yyy5a*rx(ipp,j)
2866 $ -yyy5b*rx(jm,j))*yy5
2867 y5y=(yy55*(ry(jm,j)+ry(ipp,j))-yyy5a*ry(ipp,j)
2868 $ -yyy5b*ry(jm,j))*yy5
2869 y5z=(yy55*(rz(jm,j)+rz(ipp,j))-yyy5a*rz(ipp,j)
2870 $ -yyy5b*rz(jm,j))*yy5
2883 shefx(j,11)=shefx(j,11)-sx*vbetap(i,jmm)
2884 $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
2885 shefy(j,11)=shefy(j,11)-sy*vbetap(i,jmm)
2886 $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
2887 shefz(j,11)=shefz(j,11)-sz*vbetap(i,jmm)
2888 $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
2890 ! shefx(j,11)=shefx(j,11)
2891 ! $ -sx1*vbetap1(i,jmm)-sx2*vbetap2(i,jmm)
2892 ! shefy(j,11)=shefy(j,11)
2893 ! $ -sy1*vbetap1(i,jmm)-sy2*vbetap2(i,jmm)
2894 ! shefz(j,11)=shefz(j,11)
2895 ! $ -sz1*vbetap1(i,jmm)-sz2*vbetap2(i,jmm)
2897 yy6=-(dis(ip,j)-uldhb)/dldhb
2898 y6x=rx(ip,j)/dis(ip,j)
2899 y6y=ry(ip,j)/dis(ip,j)
2900 y6z=rz(ip,j)/dis(ip,j)
2905 yy88=1.0D0/(dis(ip,j)*dis(ip,ipp))
2906 yyy8=pina1(i,jmm)/(dis(ip,j)**2)
2907 yy8=-pina1(i,jmm)/dshe
2908 y8x=(yy88*rx(ip,ipp)-yyy8*rx(ip,j))*yy8
2909 y8y=(yy88*ry(ip,ipp)-yyy8*ry(ip,j))*yy8
2910 y8z=(yy88*rz(ip,ipp)-yyy8*rz(ip,j))*yy8
2912 yy99=1.0D0/(dis(ip,j)*dis(jm,j))
2913 yyy9a=pina2(i,jmm)/(dis(ip,j)**2)
2914 yyy9b=pina2(i,jmm)/(dis(jm,j)**2)
2915 yy9=-pina2(i,jmm)/dshe
2916 y9x=(yy99*(rx(jm,j)+rx(ip,j))-yyy9a*rx(ip,j)
2917 $ -yyy9b*rx(jm,j))*yy9
2918 y9y=(yy99*(ry(jm,j)+ry(ip,j))-yyy9a*ry(ip,j)
2919 $ -yyy9b*ry(jm,j))*yy9
2920 y9z=(yy99*(rz(jm,j)+rz(ip,j))-yyy9a*rz(ip,j)
2921 $ -yyy9b*rz(jm,j))*yy9
2923 yy1010=1.0D0/(dis(jm,ipp)*dis(jm,j))
2924 yyy10=pina4(i,jmm)/(dis(jm,j)**2)
2925 yy10=-pina4(i,jmm)/dshe
2926 y10x=(yy1010*rx(jm,ipp)-yyy10*rx(jm,j))*yy10
2927 y10y=(yy1010*ry(jm,ipp)-yyy10*ry(jm,j))*yy10
2928 y10z=(yy1010*rz(jm,ipp)-yyy10*rz(jm,j))*yy10
2930 sx=y66x+y8x+y9x+y10x
2931 sy=y66y+y8y+y9y+y10y
2932 sz=y66z+y8z+y9z+y10z
2941 shefx(j,11)=shefx(j,11)-sx*vbetam(i,jmm)
2942 $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
2943 shefy(j,11)=shefy(j,11)-sy*vbetam(i,jmm)
2944 $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
2945 shefz(j,11)=shefz(j,11)-sz*vbetam(i,jmm)
2946 $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
2948 ! shefx(j,11)=shefx(j,11)
2949 ! $ -sx1*vbetam1(i,jmm)-sx2*vbetam2(i,jmm)
2950 ! shefy(j,11)=shefy(j,11)
2951 ! $ -sy1*vbetam1(i,jmm)-sy2*vbetam2(i,jmm)
2952 ! shefz(j,11)=shefz(j,11)
2953 ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
2960 c-----------------------------------------------------------------------
2961 subroutine sheetforce12
2964 parameter(maxca=800)
2965 cc**********************************************************************
2966 real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
2967 real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
2968 real*8 vbetap2(maxca,maxca),vbetam2(maxca,maxca)
2969 real*8 pin1(maxca,maxca),pin2(maxca,maxca)
2970 real*8 pin3(maxca,maxca),pin4(maxca,maxca)
2971 real*8 pina1(maxca,maxca),pina2(maxca,maxca)
2972 real*8 pina3(maxca,maxca),pina4(maxca,maxca)
2973 real*8 rx(maxca,maxca)
2974 real*8 ry(maxca,maxca),rz(maxca,maxca)
2975 real*8 bx(maxca),by(maxca),bz(maxca)
2976 real*8 dis(maxca,maxca)
2977 real*8 shefx(maxca,12),shefy(maxca,12)
2978 real*8 shefz(maxca,12)
2979 real*8 dp45,dm45,w_beta
2980 real*8 dca,dlhb,ulhb,dshe,dldhb,uldhb,
2981 $ c00,s00,ulnex,dnex
2982 integer inb,nmax,iselect
2983 cc**********************************************************************
2984 common /phys1/ inb,nmax,iselect
2986 common /difvec/ rx,ry,rz
2987 common /sheparm/ dca,dlhb,ulhb,dshe,dldhb,uldhb,
2988 $ c00,s00,ulnex,dnex
2989 common /sheconst/ dp45,dm45,w_beta
2990 common /she/ vbetap,vbetam,vbetap1,vbetap2,vbetam1,vbetam2
2991 common /shepin/ pin1,pin2,pin3,pin4,pina1,pina2,pina3,pina4
2992 common /shef/ shefx,shefy,shefz
2993 cc**********************************************************************
2995 integer j,jm,jmm,ip,i,ipp,jp
2996 real*8 yy1,y1x,y1y,y1z,y11x,y11y,y11z,yy33,yyy3,yy3,y3x,y3y,y3z
2997 real*8 yy44,yyy4a,yyy4b,yy4,y4x,y4y,y4z,yy55,yyy5,yy5,y5x,y5y,y5z
2998 real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
2999 real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
3000 real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
3001 !c*************************************************************************c
3009 yy1=-(dis(ip,j)-ulhb)/dlhb
3010 y1x=rx(ip,j)/dis(ip,j)
3011 y1y=ry(ip,j)/dis(ip,j)
3012 y1z=rz(ip,j)/dis(ip,j)
3017 yy33=1.0D0/(dis(ip,j)*dis(ip,ipp))
3018 yyy3=pin1(i,jm)/(dis(ip,j)**2)
3019 yy3=-pin1(i,jm)/dshe
3020 y3x=(yy33*rx(ip,ipp)-yyy3*rx(ip,j))*yy3
3021 y3y=(yy33*ry(ip,ipp)-yyy3*ry(ip,j))*yy3
3022 y3z=(yy33*rz(ip,ipp)-yyy3*rz(ip,j))*yy3
3023 yy44=1.0D0/(dis(ip,j)*dis(j,jp))
3025 yyy4a=pin2(i,jm)/(dis(ip,j)**2)
3026 yyy4b=pin2(i,jm)/(dis(j,jp)**2)
3027 yy4=-pin2(i,jm)/dshe
3028 y4x=(yy44*(rx(j,jp)-rx(ip,j))-yyy4a*rx(ip,j)
3029 $ +yyy4b*rx(j,jp))*yy4
3030 y4y=(yy44*(ry(j,jp)-ry(ip,j))-yyy4a*ry(ip,j)
3031 $ +yyy4b*ry(j,jp))*yy4
3032 y4z=(yy44*(rz(j,jp)-rz(ip,j))-yyy4a*rz(ip,j)
3033 $ +yyy4b*rz(j,jp))*yy4
3035 yy55=1.0D0/(dis(ipp,jp)*dis(j,jp))
3036 yyy5=pin4(i,jm)/(dis(j,jp)**2)
3037 yy5=-pin4(i,jm)/dshe
3038 y5x=(-yy55*rx(ipp,jp)+yyy5*rx(j,jp))*yy5
3039 y5y=(-yy55*ry(ipp,jp)+yyy5*ry(j,jp))*yy5
3040 y5z=(-yy55*rz(ipp,jp)+yyy5*rz(j,jp))*yy5
3053 shefx(j,12)=shefx(j,12)-sx*vbetap(i,jm)
3054 $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3055 shefy(j,12)=shefy(j,12)-sy*vbetap(i,jm)
3056 $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3057 shefz(j,12)=shefz(j,12)-sz*vbetap(i,jm)
3058 $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3060 ! shefx(j,12)=shefx(j,12)
3061 ! $ -sx1*vbetap1(i,jm)-sx2*vbetap2(i,jm)
3062 ! shefy(j,12)=shefy(j,12)
3063 ! $ -sy1*vbetap1(i,jm)-sy2*vbetap2(i,jm)
3064 ! shefz(j,12)=shefz(j,12)
3065 ! $ -sz1*vbetap1(i,jm)-sz2*vbetap2(i,jm)
3067 yy6=-(dis(ipp,j)-uldhb)/dldhb
3068 y6x=rx(ipp,j)/dis(ipp,j)
3069 y6y=ry(ipp,j)/dis(ipp,j)
3070 y6z=rz(ipp,j)/dis(ipp,j)
3075 yy88=1.0D0/(dis(ip,jp)*dis(j,jp))
3076 yyy8=pina2(i,jm)/(dis(j,jp)**2)
3077 yy8=-pina2(i,jm)/dshe
3078 y8x=(-yy88*rx(ip,jp)+yyy8*rx(j,jp))*yy8
3079 y8y=(-yy88*ry(ip,jp)+yyy8*ry(j,jp))*yy8
3080 y8z=(-yy88*rz(ip,jp)+yyy8*rz(j,jp))*yy8
3082 yy99=1.0D0/(dis(j,ipp)*dis(ip,ipp))
3083 yyy9=pina3(i,jm)/(dis(j,ipp)**2)
3084 yy9=-pina3(i,jm)/dshe
3085 y9x=(-yy99*rx(ip,ipp)+yyy9*rx(j,ipp))*yy9
3086 y9y=(-yy99*ry(ip,ipp)+yyy9*ry(j,ipp))*yy9
3087 y9z=(-yy99*rz(ip,ipp)+yyy9*rz(j,ipp))*yy9
3089 yy1010=1.0D0/(dis(j,ipp)*dis(j,jp))
3090 yyy10a=pina4(i,jm)/(dis(j,ipp)**2)
3091 yyy10b=pina4(i,jm)/(dis(j,jp)**2)
3092 yy10=-pina4(i,jm)/dshe
3093 y10x=(-yy1010*(rx(j,ipp)+rx(j,jp))+yyy10a*rx(j,ipp)
3094 $ +yyy10b*rx(j,jp))*yy10
3095 y10y=(-yy1010*(ry(j,ipp)+ry(j,jp))+yyy10a*ry(j,ipp)
3096 $ +yyy10b*ry(j,jp))*yy10
3097 y10z=(-yy1010*(rz(j,ipp)+rz(j,jp))+yyy10a*rz(j,ipp)
3098 $ +yyy10b*rz(j,jp))*yy10
3100 sx=y66x+y8x+y9x+y10x
3101 sy=y66y+y8y+y9y+y10y
3102 sz=y66z+y8z+y9z+y10z
3111 shefx(j,12)=shefx(j,12)-sx*vbetam(i,jm)
3112 $ -sx1*vbetam1(i,jm)-sx2*vbetam2(i,jm)
3113 shefy(j,12)=shefy(j,12)-sy*vbetam(i,jm)
3114 $ -sy1*vbetam1(i,jm)-sy2*vbetam2(i,jm)
3115 shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
3116 $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
3123 C===============================================================================