X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=6bb718af0a69d38643551629213163eb64393cb5;hb=da7cb646f5dae1220a85ea1473affef174d32243;hp=59b09b7b31b594c9c205b6922b16ae1e6162e6ef;hpb=65b670d78bae7c2df18f3465fdad7b4f13490e33;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 59b09b7..6bb718a 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -11,15 +11,27 @@ ! implicit none !----------------------------------------------------------------------------- +! Max. number of contacts per residue +! integer :: maxconts +!----------------------------------------------------------------------------- +! Max. number of derivatives of virtual-bond and side-chain vectors in theta +! or phi. +! integer :: maxdim +!----------------------------------------------------------------------------- +! Max. number of SC contacts +! integer :: maxcont +!----------------------------------------------------------------------------- ! Max. number of variables integer :: maxvar !----------------------------------------------------------------------------- -! Max number of torsional terms in SCCOR - integer,parameter :: maxterm_sccor=6 +! Max number of torsional terms in SCCOR in control_data +! integer,parameter :: maxterm_sccor=6 !----------------------------------------------------------------------------- ! Maximum number of SC local term fitting function coefficiants integer,parameter :: maxsccoef=65 !----------------------------------------------------------------------------- +! commom.calc common/calc/ +!----------------------------------------------------------------------------- ! commom.contacts ! common /contacts/ ! Change 12/1/95 - common block CONTACTS1 included. @@ -274,29 +286,48 @@ ! ! Compute the side-chain and electrostatic interaction energy ! - goto (101,102,103,104,105,106) ipot +! goto (101,102,103,104,105,106) ipot + select case(ipot) ! Lennard-Jones potential. - 101 call elj(evdw) +! 101 call elj(evdw) + case (1) + call elj(evdw) !d print '(a)','Exit ELJcall el' - goto 107 +! goto 107 ! Lennard-Jones-Kihara potential (shifted). - 102 call eljk(evdw) - goto 107 +! 102 call eljk(evdw) + case (2) + call eljk(evdw) +! goto 107 ! Berne-Pechukas potential (dilated LJ, angular dependence). - 103 call ebp(evdw) - goto 107 +! 103 call ebp(evdw) + case (3) + call ebp(evdw) +! goto 107 ! Gay-Berne potential (shifted LJ, angular dependence). - 104 call egb(evdw) - goto 107 +! 104 call egb(evdw) + case (4) + call egb(evdw) +! goto 107 ! Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). - 105 call egbv(evdw) - goto 107 +! 105 call egbv(evdw) + case (5) + call egbv(evdw) +! goto 107 ! Soft-sphere potential - 106 call e_softsphere(evdw) +! 106 call e_softsphere(evdw) + case (6) + call e_softsphere(evdw) ! ! Calculate electrostatic (H-bonding) energy of the main chain. ! - 107 continue +! 107 continue + case default + write(iout,*)"Wrong ipot" +! return +! 50 continue + end select +! continue !mc !mc Sep-06: egb takes care of dynamic ss bonds too @@ -324,6 +355,7 @@ .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then #endif call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) +! write (iout,*) "ELEC calc" else ees=0.0d0 evdw1=0.0d0 @@ -341,6 +373,7 @@ ! Calculate excluded-volume interaction energy between peptide groups ! and side chains. ! +!elwrite(iout,*) "in etotal calc exc;luded",ipot if (ipot.lt.6) then if(wscp.gt.0d0) then @@ -353,16 +386,20 @@ ! write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +!elwrite(iout,*) "in etotal before ebond",ipot ! ! Calculate the bond-stretching energy ! call ebond(estr) +!elwrite(iout,*) "in etotal afer ebond",ipot + ! ! Calculate the disulfide-bridge and other energy and the contributions ! from other distance constraints. - print *,'Calling EHPB' +! print *,'Calling EHPB' call edis(ehpb) +!elwrite(iout,*) "in etotal afer edis",ipot ! print *,'EHPB exitted succesfully.' ! ! Calculate the virtual-bond-angle energy. @@ -377,6 +414,7 @@ ! Calculate the SC local energy. ! call esc(escloc) +!elwrite(iout,*) "in etotal afer esc",ipot ! print *,"Processor",myrank," computed USC" ! ! Calculate the virtual-bond torsional energy. @@ -392,6 +430,7 @@ ! ! 6/23/01 Calculate double-torsional energy ! +!elwrite(iout,*) "in etotal",ipot if (wtor_d.gt.0) then call etor_d(etors_d) else @@ -423,22 +462,28 @@ ecorr6=0.0d0 eturn6=0.0d0 endif +!elwrite(iout,*) "in etotal",ipot if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1) !d write (iout,*) "multibody_hb ecorr",ecorr endif +!elwrite(iout,*) "afeter multibody hb" ! print *,"Processor",myrank," computed Ucorr" ! ! If performing constraint dynamics, call the constraint energy ! after the equilibration time if(usampl.and.totT.gt.eq_time) then +!elwrite(iout,*) "afeter multibody hb" call EconstrQ +!elwrite(iout,*) "afeter multibody hb" call Econstr_back +!elwrite(iout,*) "afeter multibody hb" else Uconst=0.0d0 Uconst_back=0.0d0 endif +!elwrite(iout,*) "after Econstr" #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 @@ -490,6 +535,8 @@ #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 #endif +!el call enerprint(energia) +!elwrite(iout,*)"finish etotal" return end subroutine etotal !----------------------------------------------------------------------------- @@ -525,7 +572,7 @@ integer :: ierr real(kind=8) :: time00 if (nfgtasks.gt.1 .and. reduce) then -!el #define DEBUG + #ifdef DEBUG write (iout,*) "energies before REDUCE" call enerprint(energia) @@ -614,7 +661,7 @@ #ifdef MPI endif #endif -!el #undef DUBUG +! call enerprint(energia) call flush(iout) return end subroutine sum_energy @@ -631,33 +678,56 @@ real(kind=8) :: kfac=2.4d0 real(kind=8) :: x,x2,x3,x4,x5,licznik=1.12692801104297249644 !el local variables - real(kind=8) :: t_bath,facT,facT2,facT3,facT4,facT5 + real(kind=8) :: t_bath,facT(6) !,facT2,facT3,facT4,facT5,facT6 + real(kind=8) :: T0=3.0d2 integer :: ierror ! facT=temp0/t_bath ! facT=2*temp0/(t_bath+temp0) if (rescale_mode.eq.0) then - facT=1.0d0 - facT2=1.0d0 - facT3=1.0d0 - facT4=1.0d0 - facT5=1.0d0 + facT(1)=1.0d0 + facT(2)=1.0d0 + facT(3)=1.0d0 + facT(4)=1.0d0 + facT(5)=1.0d0 + facT(6)=1.0d0 else if (rescale_mode.eq.1) then - facT=kfac/(kfac-1.0d0+t_bath/temp0) - facT2=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2) - facT3=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3) - facT4=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4) - facT5=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5) + facT(1)=kfac/(kfac-1.0d0+t_bath/temp0) + facT(2)=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2) + facT(3)=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3) + facT(4)=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4) + facT(5)=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5) +#ifdef WHAM_RUN +!#if defined(WHAM_RUN) || defined(CLUSTER) +#if defined(FUNCTH) +! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3) + facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0 +#elif defined(FUNCT) + facT(6)=t_bath/T0 +#else + facT(6)=1.0d0 +#endif +#endif else if (rescale_mode.eq.2) then x=t_bath/temp0 x2=x*x x3=x2*x x4=x3*x x5=x4*x - facT=licznik/dlog(dexp(x)+dexp(-x)) - facT2=licznik/dlog(dexp(x2)+dexp(-x2)) - facT3=licznik/dlog(dexp(x3)+dexp(-x3)) - facT4=licznik/dlog(dexp(x4)+dexp(-x4)) - facT5=licznik/dlog(dexp(x5)+dexp(-x5)) + facT(1)=licznik/dlog(dexp(x)+dexp(-x)) + facT(2)=licznik/dlog(dexp(x2)+dexp(-x2)) + facT(3)=licznik/dlog(dexp(x3)+dexp(-x3)) + facT(4)=licznik/dlog(dexp(x4)+dexp(-x4)) + facT(5)=licznik/dlog(dexp(x5)+dexp(-x5)) +#ifdef WHAM_RUN +!#if defined(WHAM_RUN) || defined(CLUSTER) +#if defined(FUNCTH) + facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0 +#elif defined(FUNCT) + facT(6)=t_bath/T0 +#else + facT(6)=1.0d0 +#endif +#endif else write (iout,*) "Wrong RESCALE_MODE",rescale_mode write (*,*) "Wrong RESCALE_MODE",rescale_mode @@ -666,17 +736,17 @@ #endif stop 555 endif - welec=weights(3)*fact - wcorr=weights(4)*fact3 - wcorr5=weights(5)*fact4 - wcorr6=weights(6)*fact5 - wel_loc=weights(7)*fact2 - wturn3=weights(8)*fact2 - wturn4=weights(9)*fact3 - wturn6=weights(10)*fact5 - wtor=weights(13)*fact - wtor_d=weights(14)*fact2 - wsccor=weights(21)*fact + welec=weights(3)*fact(1) + wcorr=weights(4)*fact(3) + wcorr5=weights(5)*fact(4) + wcorr6=weights(6)*fact(5) + wel_loc=weights(7)*fact(2) + wturn3=weights(8)*fact(2) + wturn4=weights(9)*fact(3) + wturn6=weights(10)*fact(5) + wtor=weights(13)*fact(1) + wtor_d=weights(14)*fact(2) + wsccor=weights(21)*fact(1) return end subroutine rescale_weights @@ -1202,7 +1272,7 @@ integer :: iint,itypi,itypi1,itypj real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi real(kind=8) :: evdw,sig0ij - + integer :: ii !cccc energy_dec=.false. ! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 @@ -1233,14 +1303,16 @@ evdw=evdw+evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,' ss' ELSE !el ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) -! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, -! & 1.0d0/vbld(j+nres) +! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,& +! 1.0d0/vbld(j+nres) !d ! write (iout,*) "i",i," j", j," itype",itype(i),itype(j) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) @@ -1269,8 +1341,12 @@ dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) ! write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi -! write (iout,*) "j",j," dc_norm", -! & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +! write (iout,*) "j",j," dc_norm",& !d +! dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +! write(iout,*)"rrij ",rrij +! write(iout,*)"xj yj zj ", xj, yj, zj +! write(iout,*)"xi yi zi ", xi, yi, zi +! write(iout,*)"c ", c(1,:), c(2,:), c(3,:) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) ! Calculate angle-dependent terms of energy and contributions to their @@ -1279,6 +1355,8 @@ sigsq=1.0D0/sigsq sig=sig0ij*dsqrt(sigsq) rij_shift=1.0D0/rij-sig+sig0ij +! write(iout,*)" rij_shift",rij_shift," rij",rij," sig",sig,& +! "sig0ij",sig0ij ! for diagnostics; uncomment ! rij_shift=1.2*sig0ij ! I hate to put IF's in the loops, but here don't have another choice!!!! @@ -1298,8 +1376,9 @@ evdwij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt -! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& -! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +! write(iout,*)"aa, bb ",aa(:,:),bb(:,:) +! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d +! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d evdwij=evdwij*eps2rt*eps3rt evdw=evdw+evdwij if (lprn) then @@ -1314,7 +1393,9 @@ endif if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & - 'evdw',i,j,evdwij + 'evdw',i,j,evdwij !,"egb" +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij ! Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1952,56 +2033,11 @@ integer :: i,iti1,iti,k,l real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2 -! allocate(Ug(2,2,nres)) !(2,2,maxres) -! allocate(Ug2(2,2,nres)) !(2,2,maxres) -! allocate(Ugder(2,2,nres)) !(2,2,maxres) -! allocate(Ug2der(2,2,nres)) !(2,2,maxres) -! allocate(obrot(2,nres)) !(2,maxres) -! allocate(obrot2(2,nres)) !(2,maxres) -! allocate(obrot_der(2,nres)) !(2,maxres) -! allocate(obrot2_der(2,nres)) !(2,maxres) -! allocate(costab2(nres)) !(maxres) -! allocate(sintab2(nres)) !(maxres) -! allocate(costab(nres)) !(maxres) -! allocate(sintab(nres)) !(maxres) - -! allocate(Ub2(2,nres)) !(2,maxres) -! allocate(Ctobr(2,nres)) !(2,maxres) -! allocate(Dtobr2(2,nres)) !(2,maxres) -! allocate(mu(2,nres)) !(2,maxres) -! allocate(muder(2,nres)) !(2,maxres) -! allocate(Ub2der(2,nres)) !(2,maxres) -! allocate(Ctobrder(2,nres)) !(2,maxres) -! allocate(Dtobr2der(2,nres)) !(2,maxres) - -! allocate(EUg(2,2,nres)) !(2,2,maxres) -! allocate(CUg(2,2,nres)) !(2,2,maxres) -! allocate(DUg(2,2,nres)) !(2,2,maxres) -! allocate(DtUg2(2,2,nres)) !(2,2,maxres) -! allocate(EUgder(2,2,nres)) !(2,2,maxres) -! allocate(CUgder(2,2,nres)) !(2,2,maxres) -! allocate(DUgder(2,2,nres)) !(2,2,maxres) -! allocate(Dtug2der(2,2,nres)) !(2,2,maxres) - -! allocate(Ug2Db1t(2,nres)) !(2,maxres) -! allocate(Ug2Db1tder(2,nres)) !(2,maxres) -! allocate(CUgb2(2,nres)) !(2,maxres) -! allocate(CUgb2der(2,nres)) !(2,maxres) - -! allocate(EUgC(2,2,nres)) !(2,2,maxres) -! allocate(EUgCder(2,2,nres)) !(2,2,maxres) -! allocate(EUgD(2,2,nres)) !(2,2,maxres) -! allocate(EUgDder(2,2,nres)) !(2,2,maxres) -! allocate(DtUg2EUg(2,2,nres)) !(2,2,maxres) -! allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres) - -! allocate(Ug2DtEUgder(2,2,2,nres)) !(2,2,2,maxres) -! allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres) - ! ! Compute the virtual-bond-torsional-angle dependent quantities needed ! to calculate the el-loc multibody terms of various order. ! +!AL el mu=0.0d0 #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else @@ -2132,7 +2168,9 @@ do k=1,2 mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) enddo -!d write (iout,*) 'mu ',mu(:,i-2) +! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2) +! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1) +! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2) !d write (iout,*) 'mu1',mu1(:,i-2) !d write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) & @@ -2171,7 +2209,6 @@ enddo endif #if defined(MPI) && defined(PARMAT) -!el #define DUBUG #ifdef DEBUG ! if (fg_rank.eq.0) then write (iout,*) "Arrays UG and UGDER before GATHER" @@ -2428,7 +2465,6 @@ !d enddo !d enddo return -!el #undef DUBUG end subroutine set_matrices !----------------------------------------------------------------------------- subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) @@ -2500,7 +2536,26 @@ !d enddo !d call check_vecgrad !d stop +! ees=0.0d0 !AS +! evdw1=0.0d0 +! eel_loc=0.0d0 +! eello_turn3=0.0d0 +! eello_turn4=0.0d0 + t_eelecij=0.0d0 + ees=0.0D0 + evdw1=0.0D0 + eel_loc=0.0d0 + eello_turn3=0.0d0 + eello_turn4=0.0d0 +! + if (icheckgrad.eq.1) then +!el +! do i=0,2*nres+2 +! dc_norm(1,i)=0.0d0 +! dc_norm(2,i)=0.0d0 +! dc_norm(3,i)=0.0d0 +! enddo do i=1,nres-1 fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i))) do k=1,3 @@ -3068,10 +3123,12 @@ ! Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij +! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33 +! if (energy_dec) write (iout,*) "muij",muij ! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) eel_loc=eel_loc+eel_loc_ij @@ -3125,6 +3182,7 @@ r0ij=2.20D0*rpp(iteli,itelj) ! r0ij=1.55D0*rpp(iteli,itelj) call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont) +!elwrite(iout,*) "num_conti",num_conti, "maxconts",maxconts if (fcont.gt.0.0D0) then num_conti=num_conti+1 if (num_conti.gt.maxconts) then @@ -4188,7 +4246,7 @@ if (i.gt.3 .and. itype(i-2).ne.ntyp1) then #ifdef OSF - phii=phi(i) + phii=phi(i) if (phii.ne.phii) phii=150.0 #else phii=phi(i) @@ -4201,7 +4259,7 @@ endif if (i.lt.nres .and. itype(i).ne.ntyp1) then #ifdef OSF - phii1=phi(i+1) + phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@ -4917,7 +4975,6 @@ sumene1x,sumene2x,sumene3x,sumene4x,& sumene1y,sumene2y,sumene3y,sumene4y,cossc,cossc1,& cosfac2xx,sinfac2yy -!el #define DEBUG #ifdef DEBUG real(kind=8) :: aincr,xxsave,sumenep,de_dxx_num,yysave,& de_dyy_num,zzsave,de_dzz_num,costsave,sintsave,& @@ -5144,7 +5201,6 @@ #ifdef DEBUG write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i) #endif -!#undef DEBUG ! ! cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres)) @@ -5219,7 +5275,6 @@ 1 continue enddo -!el #undef DUBUG return end subroutine esc !----------------------------------------------------------------------------- @@ -5375,8 +5430,8 @@ etors_ii=0.0D0 if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 & .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + itori=itortyp(itype(i-2)) + itori1=itortyp(itype(i-1)) phii=phi(i) gloci=0.0D0 ! Proline-Proline pair is a special case... @@ -5672,6 +5727,7 @@ esccor_ii=0.0D0 isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) + ! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1) phii=phi(i) do intertyp=1,3 !intertyp @@ -8995,12 +9051,12 @@ #endif #ifdef MPI include 'mpif.h' -!el#endif +#endif real(kind=8),dimension(3,nres) :: gradbufc,gradbufx,gradbufc_sum,& gloc_scbuf !(3,maxres) real(kind=8),dimension(4*nres) :: glocbuf !(4*maxres) -#endif +!#endif !el local variables integer :: i,j,k,ierror,ierr real(kind=8) :: gvdwc_norm,gvdwc_scp_norm,gelc_norm,gvdwpp_norm,& @@ -10199,17 +10255,13 @@ icall=0 call geom_to_var(nvar,x) if (.not.split_ene) then -write(iout,*) 'Calling CHECK_ECARTINT if' call etotal(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' etot=energia(0) !el call enerprint(energia) -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' call flush(iout) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' write (iout,*) "exit cartgrad" call flush(iout) icall =1 @@ -10219,7 +10271,6 @@ write(iout,*) 'Calling CHECK_ECARTINT if' do j=1,3 grad_s(j,0)=gcart(j,0) enddo -!elwrite(iout,*) 'Calling CHECK_ECARTINT if' do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) @@ -10227,7 +10278,6 @@ write(iout,*) 'Calling CHECK_ECARTINT if' enddo enddo else -write(iout,*) 'Calling CHECK_ECARTIN else.' !- split gradient check call zerograd call etotal_long(energia) @@ -10421,14 +10471,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' call zerograd aincr=1.0D-7 print '(a)','Calling CHECK_INT.' -write(iout,*) 'Calling CHECK_INT.' nf=0 nfl=0 icg=1 call geom_to_var(nvar,x) call var_to_geom(nvar,x) call chainbuild -write(iout,*) 'Calling CHECK_INT.' icall=1 print *,'ICG=',ICG call etotal(energia) @@ -10447,7 +10495,7 @@ write(iout,*) 'Calling CHECK_INT.' nfl=3 !d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) - write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp +!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp icall=1 do i=1,nvar xi=x(i) @@ -10485,7 +10533,6 @@ write(iout,*) 'Calling CHECK_INT.' i,key,ii,gg(i),gana(i),& 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr) enddo -write(iout,*) "jestesmy sobie w check eint!!" return end subroutine check_eint !----------------------------------------------------------------------------- @@ -11325,6 +11372,8 @@ write(iout,*) "jestesmy sobie w check eint!!" if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'evdw',i,j,evdwij +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,"egb_long" ! Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -11468,6 +11517,8 @@ write(iout,*) "jestesmy sobie w check eint!!" if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'evdw',i,j,evdwij +! if (energy_dec) write (iout,*) & +! 'evdw',i,j,evdwij,"egb_short" ! Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -12398,11 +12449,11 @@ write(iout,*) "jestesmy sobie w check eint!!" ! Contribution to the local-electrostatic energy coming from the i-j pair eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) & +a33*muij(4) -!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij +! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij -! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) +! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d eel_loc=eel_loc+eel_loc_ij ! Partial derivatives in virtual-bond dihedral angles gamma @@ -13708,7 +13759,9 @@ write(iout,*) "jestesmy sobie w check eint!!" call sum_gradient #ifdef TIMING #endif +!el write (iout,*) "After sum_gradient" #ifdef DEBUG +!el write (iout,*) "After sum_gradient" do i=1,nres-1 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) @@ -13729,10 +13782,12 @@ write(iout,*) "jestesmy sobie w check eint!!" gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i) enddo endif +!elwrite (iout,*) "After sum_gradient" #ifdef TIMING time01=MPI_Wtime() #endif call intcartderiv +!elwrite (iout,*) "After sum_gradient" #ifdef TIMING time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01 #endif @@ -14089,11 +14144,6 @@ write(iout,*) "jestesmy sobie w check eint!!" endif enddo !alculate derivative of Tauangle - do i=1,nres-1 - do j=1,3 - dc_norm2(j,i+nres)=-dc_norm(j,i+nres) - enddo - enddo #ifdef PARINTDER do i=itau_start,itau_end #else @@ -14112,7 +14162,10 @@ write(iout,*) "jestesmy sobie w check eint!!" cost=dcos(theta(i)) cost1=dcos(omicron(2,i-1)) cosg=dcos(tauangle(1,i)) +!elwrite(iout,*) " vecpr5",i,nres do j=1,3 +!elwrite(iout,*) " vecpreee",i,nres,j,i-2+nres +!elwrite(iout,*) " vecpr5",dc_norm2(1,1) dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres) ! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" enddo @@ -15027,7 +15080,7 @@ write(iout,*) "jestesmy sobie w check eint!!" integer :: kstart,kend,lstart,lend,idummy real(kind=8) :: delta=1.0d-7 !el local variables - integer :: i,ii,j + integer :: i,ii,j ! real(kind=8) :: ! For the backbone do i=0,nres-1 @@ -15273,6 +15326,7 @@ write(iout,*) "jestesmy sobie w check eint!!" real(kind=8) :: deps,ssx0,ljx0 !-------END TESTING CODE + eij=0.0d0 i=resi j=resj @@ -15532,31 +15586,31 @@ write(iout,*) "jestesmy sobie w check eint!!" endif if (havebond) then -#ifndef CLUST -#ifndef WHAM +!#ifndef CLUST +!#ifndef WHAM ! if (dyn_ssbond_ij(i,j).eq.1.0d300) then ! write(iout,'(a15,f12.2,f8.1,2i5)') ! & "SSBOND_E_FORM",totT,t_bath,i,j ! endif -#endif -#endif +!#endif +!#endif dyn_ssbond_ij(i,j)=eij else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then dyn_ssbond_ij(i,j)=1.0d300 -#ifndef CLUST -#ifndef WHAM +!#ifndef CLUST +!#ifndef WHAM ! write(iout,'(a15,f12.2,f8.1,2i5)') ! & "SSBOND_E_BREAK",totT,t_bath,i,j -#endif -#endif +!#endif +!#endif endif !-------TESTING CODE - if (checkstop) then +!el if (checkstop) then if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') & "CHECKSTOP",rij,eij,ed echeck(jcheck)=eij - endif +!el endif enddo if (checkstop) then write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps @@ -15649,11 +15703,7 @@ write(iout,*) "jestesmy sobie w check eint!!" ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.SETUP' -#ifndef CLUST -#ifndef WHAM ! include 'COMMON.MD' -#endif -#endif ! Local variables real(kind=8) :: emin integer :: i,j,imin,ierr @@ -15788,7 +15838,7 @@ write(iout,*) "jestesmy sobie w check eint!!" !----------------------------------------------------------------------------- #ifdef WHAM subroutine read_ssHist - implicit none +! implicit none ! Includes ! include 'DIMENSIONS' ! include "DIMENSIONS.FREE" @@ -15844,28 +15894,42 @@ write(iout,*) "jestesmy sobie w check eint!!" maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond !---------------------- ! arrays in subroutine init_int_table +!el#ifdef MPI +!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1) +!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1) +!el#endif allocate(nint_gr(nres)) allocate(nscp_gr(nres)) allocate(ielstart(nres)) - allocate(ielend(nres)) !(maxres) + allocate(ielend(nres)) +!(maxres) allocate(istart(nres,maxint_gr)) - allocate(iend(nres,maxint_gr)) !(maxres,maxint_gr) + allocate(iend(nres,maxint_gr)) +!(maxres,maxint_gr) allocate(iscpstart(nres,maxint_gr)) - allocate(iscpend(nres,maxint_gr)) !(maxres,maxint_gr) + allocate(iscpend(nres,maxint_gr)) +!(maxres,maxint_gr) allocate(ielstart_vdw(nres)) - allocate(ielend_vdw(nres)) !(maxres) + allocate(ielend_vdw(nres)) +!(maxres) - allocate(lentyp(0:nfgtasks-1)) !(0:maxprocs-1) + allocate(lentyp(0:nfgtasks-1)) +!(0:maxprocs-1) !---------------------- ! commom.contacts ! common /contacts/ if(.not.allocated(icont_ref)) allocate(icont_ref(2,maxcont)) - allocate(icont(2,maxcont)) !(2,maxcont) + allocate(icont(2,maxcont)) +!(2,maxcont) ! common /contacts1/ - allocate(num_cont(0:nres+4)) !(maxres) - allocate(jcont(maxconts,nres)) !(maxconts,maxres) - allocate(facont(maxconts,nres)) !(maxconts,maxres) - allocate(gacont(3,maxconts,nres)) !(3,maxconts,maxres) + allocate(num_cont(0:nres+4)) +!(maxres) + allocate(jcont(maxconts,nres)) +!(maxconts,maxres) + allocate(facont(maxconts,nres)) +!(maxconts,maxres) + allocate(gacont(3,maxconts,nres)) +!(3,maxconts,maxres) ! common /contacts_hb/ allocate(gacontp_hb1(3,maxconts,nres)) allocate(gacontp_hb2(3,maxconts,nres)) @@ -15874,31 +15938,42 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(gacontm_hb2(3,maxconts,nres)) allocate(gacontm_hb3(3,maxconts,nres)) allocate(gacont_hbr(3,maxconts,nres)) - allocate(grij_hb_cont(3,maxconts,nres)) !(3,maxconts,maxres) + allocate(grij_hb_cont(3,maxconts,nres)) +!(3,maxconts,maxres) allocate(facont_hb(maxconts,nres)) allocate(ees0p(maxconts,nres)) allocate(ees0m(maxconts,nres)) - allocate(d_cont(maxconts,nres)) !(maxconts,maxres) - allocate(num_cont_hb(nres)) !(maxres) - allocate(jcont_hb(maxconts,nres)) !(maxconts,maxres) + allocate(d_cont(maxconts,nres)) +!(maxconts,maxres) + allocate(num_cont_hb(nres)) +!(maxres) + allocate(jcont_hb(maxconts,nres)) +!(maxconts,maxres) ! common /rotat/ allocate(Ug(2,2,nres)) allocate(Ugder(2,2,nres)) allocate(Ug2(2,2,nres)) - allocate(Ug2der(2,2,nres)) !(2,2,maxres) + allocate(Ug2der(2,2,nres)) +!(2,2,maxres) allocate(obrot(2,nres)) allocate(obrot2(2,nres)) allocate(obrot_der(2,nres)) - allocate(obrot2_der(2,nres)) !(2,maxres) + allocate(obrot2_der(2,nres)) +!(2,maxres) ! common /precomp1/ allocate(mu(2,nres)) allocate(muder(2,nres)) allocate(Ub2(2,nres)) + do i=1,nres + Ub2(1,i)=0.0d0 + Ub2(2,i)=0.0d0 + enddo allocate(Ub2der(2,nres)) allocate(Ctobr(2,nres)) allocate(Ctobrder(2,nres)) allocate(Dtobr2(2,nres)) - allocate(Dtobr2der(2,nres)) !(2,maxres) + allocate(Dtobr2der(2,nres)) +!(2,maxres) allocate(EUg(2,2,nres)) allocate(EUgder(2,2,nres)) allocate(CUg(2,2,nres)) @@ -15906,25 +15981,30 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(DUg(2,2,nres)) allocate(Dugder(2,2,nres)) allocate(DtUg2(2,2,nres)) - allocate(DtUg2der(2,2,nres)) !(2,2,maxres) + allocate(DtUg2der(2,2,nres)) +!(2,2,maxres) ! common /precomp2/ allocate(Ug2Db1t(2,nres)) allocate(Ug2Db1tder(2,nres)) allocate(CUgb2(2,nres)) - allocate(CUgb2der(2,nres)) !(2,maxres) + allocate(CUgb2der(2,nres)) +!(2,maxres) allocate(EUgC(2,2,nres)) allocate(EUgCder(2,2,nres)) allocate(EUgD(2,2,nres)) allocate(EUgDder(2,2,nres)) allocate(DtUg2EUg(2,2,nres)) - allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres) + allocate(Ug2DtEUg(2,2,nres)) +!(2,2,maxres) allocate(Ug2DtEUgder(2,2,2,nres)) - allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres) + allocate(DtUg2EUgder(2,2,2,nres)) +!(2,2,2,maxres) ! common /rotat_old/ allocate(costab(nres)) allocate(sintab(nres)) allocate(costab2(nres)) - allocate(sintab2(nres)) !(maxres) + allocate(sintab2(nres)) +!(maxres) ! common /dipmat/ allocate(a_chuj(2,2,maxconts,nres)) !(2,2,maxconts,maxres)(maxconts=maxres/4) @@ -15934,24 +16014,33 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(ncont_sent(nres)) allocate(ncont_recv(nres)) - allocate(iat_sent(nres)) !(maxres) + allocate(iat_sent(nres)) +!(maxres) allocate(iint_sent(4,nres,nres)) - allocate(iint_sent_local(4,nres,nres)) !(4,maxres,maxres) + allocate(iint_sent_local(4,nres,nres)) +!(4,maxres,maxres) allocate(iturn3_sent(4,0:nres+4)) allocate(iturn4_sent(4,0:nres+4)) allocate(iturn3_sent_local(4,nres)) - allocate(iturn4_sent_local(4,nres)) !(4,maxres) + allocate(iturn4_sent_local(4,nres)) +!(4,maxres) allocate(itask_cont_from(0:nfgtasks-1)) - allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1) + allocate(itask_cont_to(0:nfgtasks-1)) +!(0:max_fg_procs-1) + + !---------------------- ! commom.deriv; ! common /derivat/ allocate(dcdv(6,maxdim)) - allocate(dxdv(6,maxdim)) !(6,maxdim) - allocate(dxds(6,nres)) !(6,maxres) + allocate(dxdv(6,maxdim)) +!(6,maxdim) + allocate(dxds(6,nres)) +!(6,maxres) allocate(gradx(3,nres,0:2)) - allocate(gradc(3,nres,0:2)) !(3,maxres,2) + allocate(gradc(3,nres,0:2)) +!(3,maxres,2) allocate(gvdwx(3,nres)) allocate(gvdwc(3,nres)) allocate(gelc(3,nres)) @@ -15969,32 +16058,41 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(gcorr6_turn_long(3,nres)) allocate(gradxorr(3,nres)) allocate(gradcorr5(3,nres)) - allocate(gradcorr6(3,nres)) !(3,maxres) + allocate(gradcorr6(3,nres)) +!(3,maxres) allocate(gloc(0:maxvar,0:2)) - allocate(gloc_x(0:maxvar,2)) !(maxvar,2) + allocate(gloc_x(0:maxvar,2)) +!(maxvar,2) allocate(gel_loc(3,nres)) allocate(gel_loc_long(3,nres)) allocate(gcorr3_turn(3,nres)) allocate(gcorr4_turn(3,nres)) allocate(gcorr6_turn(3,nres)) allocate(gradb(3,nres)) - allocate(gradbx(3,nres)) !(3,maxres) + allocate(gradbx(3,nres)) +!(3,maxres) allocate(gel_loc_loc(maxvar)) allocate(gel_loc_turn3(maxvar)) allocate(gel_loc_turn4(maxvar)) allocate(gel_loc_turn6(maxvar)) allocate(gcorr_loc(maxvar)) allocate(g_corr5_loc(maxvar)) - allocate(g_corr6_loc(maxvar)) !(maxvar) + allocate(g_corr6_loc(maxvar)) +!(maxvar) allocate(gsccorc(3,nres)) - allocate(gsccorx(3,nres)) !(3,maxres) - allocate(gsccor_loc(nres)) !(maxres) - allocate(dtheta(3,2,nres)) !(3,2,maxres) + allocate(gsccorx(3,nres)) +!(3,maxres) + allocate(gsccor_loc(nres)) +!(maxres) + allocate(dtheta(3,2,nres)) +!(3,2,maxres) allocate(gscloc(3,nres)) - allocate(gsclocx(3,nres)) !(3,maxres) + allocate(gsclocx(3,nres)) +!(3,maxres) allocate(dphi(3,3,nres)) allocate(dalpha(3,3,nres)) - allocate(domega(3,3,nres)) !(3,3,maxres) + allocate(domega(3,3,nres)) +!(3,3,maxres) ! common /deriv_scloc/ allocate(dXX_C1tab(3,nres)) allocate(dYY_C1tab(3,nres)) @@ -16004,10 +16102,13 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(dZZ_Ctab(3,nres)) allocate(dXX_XYZtab(3,nres)) allocate(dYY_XYZtab(3,nres)) - allocate(dZZ_XYZtab(3,nres)) !(3,maxres) + allocate(dZZ_XYZtab(3,nres)) +!(3,maxres) ! common /mpgrad/ allocate(jgrad_start(nres)) - allocate(jgrad_end(nres)) !(maxres) + allocate(jgrad_end(nres)) +!(maxres) +!---------------------- ! common /indices/ allocate(ibond_displ(0:nfgtasks-1)) @@ -16023,20 +16124,25 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(iset_displ(0:nfgtasks-1)) allocate(iset_count(0:nfgtasks-1)) allocate(iint_count(0:nfgtasks-1)) - allocate(iint_displ(0:nfgtasks-1)) !(0:max_fg_procs-1) + allocate(iint_displ(0:nfgtasks-1)) +!(0:max_fg_procs-1) !---------------------- ! common.MD ! common /mdgrad/ allocate(gcart(3,0:nres)) - allocate(gxcart(3,0:nres)) !(3,0:MAXRES) + allocate(gxcart(3,0:nres)) +!(3,0:MAXRES) allocate(gradcag(3,nres)) - allocate(gradxag(3,nres)) !(3,MAXRES) + allocate(gradxag(3,nres)) +!(3,MAXRES) ! common /back_constr/ !el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back) allocate(dutheta(nres)) - allocate(dugamma(nres)) !(maxres) + allocate(dugamma(nres)) +!(maxres) allocate(duscdiff(3,nres)) - allocate(duscdiffx(3,nres)) !(3,maxres) + allocate(duscdiffx(3,nres)) +!(3,maxres) !el i io:read_fragments ! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20) ! allocate((:,:,:),allocatable :: ifrag_back !(3,maxfrag_back,maxprocs/20) @@ -16052,7 +16158,8 @@ write(iout,*) "jestesmy sobie w check eint!!" allocate(dUdconst(3,0:nres)) allocate(dUdxconst(3,0:nres)) allocate(dqwol(3,0:nres)) - allocate(dxqwol(3,0:nres)) !(3,0:MAXRES) + allocate(dxqwol(3,0:nres)) +!(3,0:MAXRES) !---------------------- ! common.sbridge ! common /sbridge/ in io_common: read_bridge @@ -16062,7 +16169,8 @@ write(iout,*) "jestesmy sobie w check eint!!" !el integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane ! common /dyn_ssbond/ ! and side-chain vectors in theta or phi. - allocate(dyn_ssbond_ij(0:nres+4,0:nres+4)) !(maxres,maxres) + allocate(dyn_ssbond_ij(0:nres+4,0:nres+4)) +!(maxres,maxres) do i=1,nres do j=i+1,nres dyn_ssbond_ij(i,j)=1.0d300 @@ -16070,9 +16178,11 @@ write(iout,*) "jestesmy sobie w check eint!!" enddo if (nss.gt.0) then - allocate(idssb(nss),jdssb(nss)) !(maxdim) + allocate(idssb(nss),jdssb(nss)) +!(maxdim) endif - allocate(dyn_ss_mask(nres)) !(maxres) + allocate(dyn_ss_mask(nres)) +!(maxres) do i=1,nres dyn_ss_mask(i)=.false. enddo @@ -16091,59 +16201,32 @@ write(iout,*) "jestesmy sobie w check eint!!" ! allocate(vlor2sccor(maxterm_sccor,20,20)) ! allocate(vlor3sccor(maxterm_sccor,20,20)) !(maxterm_sccor,20,20) !---------------- - allocate(gloc_sc(3,0:2*nres,0:10)) !(3,0:maxres2,10)maxres2=2*maxres + allocate(gloc_sc(3,0:2*nres,0:10)) +!(3,0:maxres2,10)maxres2=2*maxres allocate(dcostau(3,3,3,2*nres)) allocate(dsintau(3,3,3,2*nres)) allocate(dtauangle(3,3,3,2*nres)) allocate(dcosomicron(3,3,3,2*nres)) - allocate(domicron(3,3,3,2*nres)) !(3,3,3,maxres2)maxres2=2*maxres -!---------------------- -! common.scrot -! Parameters of the SC rotamers (local) term -! common/scrot/ in io_conf: parmread -! allocate((:,:),allocatable :: sc_parmin !(maxsccoef,ntyp) -!---------------------- -! common.torcnstr -! common /torcnstr/ -!el in io_conf:molread -! allocate((:),allocatable :: idih_constr,idih_nconstr !(maxdih_constr) -! allocate((:),allocatable :: phi0,drange !(maxdih_constr) -!---------------------- -! common.torsion -! common/torsion/ in io_conf: parmread -! allocate((:,:,:),allocatable :: v0 !(-maxtor:maxtor,-maxtor:maxtor,2) -! allocate((:,:,:,:),allocatable :: v1,v2 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) -! allocate((:,:,:),allocatable :: vlor1 !(maxlor,-maxtor:maxtor,-maxtor:maxtor) -! allocate((:,:,:),allocatable :: vlor2,vlor3 !(maxlor,maxtor,maxtor) -! allocate((:),allocatable :: itortyp !(-ntyp1:ntyp1) -! allocate((:,:,:),allocatable :: nterm,nlor !(-maxtor:maxtor,-maxtor:maxtor,2) -! -! common /torsiond/ in io_conf: parmread -! allocate((:,:,:,:,:,:),allocatable :: v1c,v1s - !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2) -! allocate((:,:,:,:,:,:),allocatable :: v2c,v2s - !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2) -! allocate((:,:,:,:),allocatable :: ntermd_1,ntermd_2 - !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2) -! common/fourier/ in io_conf: parmread -! allocate((:,:),allocatable :: b1,b2,& -! b1tilde !(2,-maxtor:maxtor) -! allocate((:,:,:),allocatable :: cc,dd,ee,& -! ctilde,dtilde !(2,2,-maxtor:maxtor) + allocate(domicron(3,3,3,2*nres)) +!(3,3,3,maxres2)maxres2=2*maxres !---------------------- ! common.var ! common /restr/ - allocate(varall(maxvar)) !(maxvar)(maxvar=6*maxres) + allocate(varall(maxvar)) +!(maxvar)(maxvar=6*maxres) allocate(mask_theta(nres)) allocate(mask_phi(nres)) - allocate(mask_side(nres)) !(maxres) + allocate(mask_side(nres)) +!(maxres) !---------------------- ! common.vectors ! common /vectors/ allocate(uy(3,nres)) - allocate(uz(3,nres)) !(3,maxres) + allocate(uz(3,nres)) +!(3,maxres) allocate(uygrad(3,3,2,nres)) - allocate(uzgrad(3,3,2,nres)) !(3,3,2,maxres) + allocate(uzgrad(3,3,2,nres)) +!(3,3,2,maxres) return end subroutine alloc_ener_arrays