X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fdfa.F;h=412943ac06518dc4e3fa3a20c196eb8624729515;hb=5b984996905504c1fdc0d3e94b036806f9126f7a;hp=576910c65314275a2b065d52319a9da0c30871b5;hpb=b847bacd7f528a68d640b3d58d4f5b9ce3c6ede2;p=unres.git diff --git a/source/unres/src_MD/dfa.F b/source/unres/src_MD/dfa.F index 576910c..412943a 100644 --- a/source/unres/src_MD/dfa.F +++ b/source/unres/src_MD/dfa.F @@ -74,6 +74,9 @@ C include 'COMMON.IOUNITS' include 'COMMON.CHAIN' include 'COMMON.DFA' + include 'COMMON.FFIELD' + include 'COMMON.CONTROL' + include 'COMMON.SETUP' C NOTE THAT FILENAMES are FIXED, CURRENTLY!! @@ -93,6 +96,7 @@ C open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33) goto 34 33 write(iout,'(a)') 'Error opening dist_dfa.dat file' + print *,'Error opening dist_dfa.dat file' stop 34 continue write(iout,'(a)') 'dist_dfa.dat is opened!' @@ -125,7 +129,8 @@ C READ ANGLE RESTRAINTS C PHI RESTRAINTS open(iodfa, file='phi_dfa.dat',status='old',err=35) goto 36 - 35 write(iout,'(a)') 'Error opening dist_dfa.dat file' + 35 write(iout,'(a)') 'Error opening phi_dfa.dat file' + print *, 'Error opening phi_dfa.dat file' stop 36 continue @@ -164,7 +169,8 @@ C READ NUMBER OF RESTRAINTS C THETA RESTRAINTS open(iodfa, file='theta_dfa.dat',status='old',err=41) goto 42 - 41 write(iout,'(a)') 'Error opening dist_dfa.dat file' + 41 write(iout,'(a)') 'Error opening theta_dfa.dat file' + print *,'Error opening theta_dfa.dat file' stop 42 continue write(iout,'(a)') 'theta_dfa.dat is opened!' @@ -204,6 +210,7 @@ C NUMBER OF NEIGHBOR CAs open(iodfa,file='nei_dfa.dat',status='old',err=37) goto 38 37 write(iout,'(a)') 'Error opening nei_dfa.dat file' + print *,'Error opening nei_dfa.dat file' stop 38 continue write(iout,'(a)') 'nei_dfa.dat is opened!' @@ -235,10 +242,42 @@ C write(*,*) 'READ NEI:',i,j,fnei(i,j) close(iodfa) C END OF NEIGHBORING CA +C init parallel +C BETA is not parallel ! +#ifdef MPI + if (wdfa_beta.ne.0.0 .and. nfgtasks.gt.1) then + write (iout,*) "ERRROR dfa_beta works only for FGPROCS=1" + print *,"ERRROR dfa_beta works only for FGPROCS=1" + stop + endif + call int_bounds(idfadis,idfadis_start,idfadis_end) + call int_bounds(idfaphi,idfaphi_start,idfaphi_end) + call int_bounds(idfathe,idfathe_start,idfathe_end) + call int_bounds(idfanei,idfanei_start,idfanei_end) + if (me.eq.king .or. .not. out1file) + & write (iout,*) "DFA MPI ", + & "idfadis ",idfadis,idfadis_start,idfadis_end, + & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end, + & "idfathe ",idfathe,idfathe_start,idfathe_end, + & "idfanei ",idfanei,idfanei_start,idfanei_end +#else + idfadis_start=1 + idfadis_end=idfadis + idfaphi_start=1 + idfaphi_end=idfaphi + idfathe_start=1 + idfathe_end=idfathe + idfanei_start=1 + idfanei_end=idfanei +#endif + + C READ BETA RESTRAINT + if (wdfa_beta.eq.0.0) return open(iodfa, file='beta_dfa.dat',status='old',err=39) goto 40 39 write(iout,'(a)') 'Error opening beta_dfa.dat file' + print *,'Error opening beta_dfa.dat file' stop 40 continue write(iout,'(a)') 'beta_dfa.dat is opened!' @@ -258,6 +297,7 @@ c write(*,*) 'BETA:',i,j,wtmp,wshet(i,j) close(iodfa) C END OF BETA RESTRAINT + return END @@ -278,7 +318,7 @@ C END OF BETA RESTRAINT edfadis=0 gdfad=0.0d0 - do i=1, idfadis + do i=idfadis_start,idfadis_end iatm1=idislis(1,i)+ishiftca iatm2=idislis(2,i)+ishiftca @@ -361,7 +401,7 @@ C DFA torsion angle gdfat(:,:) = 0.0d0 C START OF PHI ANGLE - do i=1, idfaphi + do i=idfaphi_start,idfaphi_end aphi = 0.0d0 do iii=1,5 @@ -560,7 +600,7 @@ c end of single assignment statement C END OF PHI RESTRAINT C START OF THETA ANGLE - do i=1, idfathe + do i=idfathe_start,idfathe_end athe = 0.0d0 do iii=1,5 @@ -794,7 +834,7 @@ C DFA neighboring CA restraint c print*, 's1:', s1(:) c print*, 's2:', s2(:) - do i=1, idfanei + do i=idfanei_start,idfanei_end kshnum=kshell(i) n1atom=ineilis(i)+ishiftca @@ -830,9 +870,9 @@ c write(*,*) n1atom, j, dist c write(*,*) 'case1:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -876,9 +916,9 @@ c neighbor atoms tmp_n=tmp_n+1.0d0 c write(*,*) 'case4:',tmp_n - t1dx=t1dx+0.0d0 - t1dy=t1dy+0.0d0 - t1dz=t1dz+0.0d0 +cc t1dx=t1dx+0.0d0 +cc t1dy=t1dy+0.0d0 +cc t1dz=t1dz+0.0d0 t2dx(j)=0.0d0 t2dy(j)=0.0d0 t2dz(j)=0.0d0 @@ -1035,8 +1075,8 @@ C End of sheet variables double precision enebet enebet=0.0d0 - bx=0.0d0;by=0.0d0;bz=0.0d0 - shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 +c bx=0.0d0;by=0.0d0;bz=0.0d0 +c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0 gdfab=0.0d0 @@ -2119,6 +2159,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc****************************************************************************** c real*8 dfaexp(15001) @@ -2152,6 +2195,9 @@ cc********************************************************************** real*8 uum, uup real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2 + real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca) + common /shetf/ shetfx,shetfy,shetfz + common /sheca/ bx,by,bz common /phys1/ inb,nmax,iselect common /kyori2/ dis @@ -2170,6 +2216,7 @@ cc********************************************************************** common /sinsin/ sth real*8 r_pair_mat(maxca,maxca) + real*8 e_gcont,fprim_gcont,de_gcont ci integer istrand(maxca,maxca) ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca) ci common /shetest/ istrand,istrand_p,istrand_m @@ -2190,6 +2237,12 @@ c do i=1,inb-7 do j=i+4,inb-3 + + if (dis(i,j).lt.dfa_cutoff) then + call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + + ip=i+1 ipp=i+2 jp=j+1 @@ -2404,7 +2457,30 @@ c uum=vbetam1(i,j)+vbetam2(i,j) vbet(i,j)=uup+uum vbetp=vbetp+uup vbetm=vbetm+uum - vbeta=vbeta+vbet(i,j) + vbeta=vbeta+vbet(i,j)*e_gcont + + + if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then +c gradient correction from gcont + de_gcont=vbet(i,j)*fprim_gcont/dis(i,j) + shetfx(i)=shetfx(i) + de_gcont*rx(i,j) + shetfy(i)=shetfy(i) + de_gcont*ry(i,j) + shetfz(i)=shetfz(i) + de_gcont*rz(i,j) + + shetfx(j)=shetfx(j) - de_gcont*rx(i,j) + shetfy(j)=shetfy(j) - de_gcont*ry(i,j) + shetfz(j)=shetfz(j) - de_gcont*rz(i,j) + +c energy correction from gcont + vbet(i,j)=vbet(i,j)*e_gcont + vbetap(i,j)=vbetap(i,j)*e_gcont + vbetap1(i,j)=vbetap1(i,j)*e_gcont + vbetap2(i,j)=vbetap2(i,j)*e_gcont + vbetam(i,j)=vbetam(i,j)*e_gcont + vbetam1(i,j)=vbetam1(i,j)*e_gcont + vbetam2(i,j)=vbetam2(i,j)*e_gcont + endif + ci elseif(istrand(i,j).eq.0)then ci vbet(i,j)=0 @@ -2423,7 +2499,16 @@ c write(*,*) 'vbetp:',vbetp c write(*,*) 'vbetm:',vbetm c write(*,*) 'vbet(i,j):',vbet(i,j) c stop - + + else + vbetap(i,j)=0 + vbetap1(i,j)=0 + vbetap2(i,j)=0 + vbetam(i,j)=0 + vbetam1(i,j)=0 + vbetam2(i,j)=0 + vbet(i,j)=0 + endif enddo enddo @@ -2447,6 +2532,9 @@ c------------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbet(maxca,maxca) real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) @@ -2505,6 +2593,7 @@ c local variables real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8 real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2 real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2 + real*8 e_gcont,fprim_gcont C-------------------------------------------------------------------------------- do i=4,inb-4 im3=i-3 @@ -2744,6 +2833,9 @@ c---------------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2784,11 +2876,17 @@ c local variables real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b + real*8 e_gcont,fprim_gcont c******************************************************************************** do i=3,inb-5 imm=i-2 im=i-1 do j=i+2,inb-3 + + if (dis(imm,j).lt.dfa_cutoff) then + call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -2912,7 +3010,8 @@ ci & .and.(istrand_p(imm,j)+istrand_m(imm,j)).ge.1) then ! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j) ! shefz(i,5)=shefz(i,5) ! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j) - + + endif ci endif enddo @@ -2925,6 +3024,9 @@ c--------------------------------------------------------------------------c implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -2965,11 +3067,17 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4 real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b + real*8 e_gcont,fprim_gcont C******************************************************************************** do i=2,inb-6 ip=i+1 im=i-1 do j=i+3,inb-3 + + if (dis(im,j).lt.dfa_cutoff) then + call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + jp=j+1 jpp=j+2 @@ -3091,7 +3199,8 @@ ci & .and.(istrand_p(im,j)+istrand_m(im,j)).ge.1) then ! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j) ! shefz(i,6)=shefz(i,6) ! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j) - + + endif ci endif enddo @@ -3104,6 +3213,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3144,12 +3256,18 @@ C local variables real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6 real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8 + real*8 e_gcont,fprim_gcont C******************************************************************************** do j=7,inb-1 jm=j-1 jmm=j-2 do i=1,j-6 + + if (dis(i,jmm).lt.dfa_cutoff) then + call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3273,6 +3391,7 @@ ci & .and.(istrand_p(i,jmm)+istrand_m(i,jmm)).ge.1) then ! shefz(j,11)=shefz(j,11) ! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm) + endif ci endif enddo @@ -3285,6 +3404,9 @@ c----------------------------------------------------------------------- implicit none integer maxca parameter(maxca=800) + real*8 dfa_cutoff,dfa_cutoff_delta + parameter(dfa_cutoff=15.5d0) + parameter(dfa_cutoff_delta=0.5d0) cc********************************************************************** real*8 vbetap(maxca,maxca),vbetam(maxca,maxca) real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca) @@ -3324,11 +3446,17 @@ C local variables real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8 + real*8 e_gcont,fprim_gcont !c*************************************************************************c do j=6,inb-2 jp=j+1 jm=j-1 do i=1,j-5 + + if (dis(i,jm).lt.dfa_cutoff) then + call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0, + & dfa_cutoff_delta,e_gcont,fprim_gcont) + ip=i+1 ipp=i+2 @@ -3445,6 +3573,8 @@ ci & .and.(istrand_p(i,jm)+istrand_m(i,jm)).ge.1) then shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm) $ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm) + endif + ci endif ENDDO