From: Adam Kazimierz Sieradzan Date: Wed, 21 May 2014 11:43:27 +0000 (-0400) Subject: diagnoza dlaczego WHAM nie dziala X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=a1bbedcfb50a59f5ae082fe9016af63119b205cc;hp=0d0c28ca14abcf99ccb778cdba15a9646736571d;p=unres.git diagnoza dlaczego WHAM nie dziala --- diff --git a/bin/wham_multparm-D-sccor-oldparm b/bin/wham_multparm-D-sccor-oldparm deleted file mode 100755 index 008abf2..0000000 Binary files a/bin/wham_multparm-D-sccor-oldparm and /dev/null differ diff --git a/source/unres/src_MD-M/MREMD.F b/source/unres/src_MD-M/MREMD.F index d55a95b..d9ddba2 100644 --- a/source/unres/src_MD-M/MREMD.F +++ b/source/unres/src_MD-M/MREMD.F @@ -525,7 +525,8 @@ c Variable time step algorithm. ugamma_cache(i,ntwx_cache)=ugamma(i) uscdiff_cache(i,ntwx_cache)=uscdiff(i) enddo - +C print *,'przed returnbox' + call returnbox do i=1,nres*2 do j=1,3 c_cache(j,i,ntwx_cache)=c(j,i) @@ -1861,4 +1862,228 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) if(me.eq.king) close(irest2) return end - +c------------------------------------------ + subroutine returnbox + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + include 'COMMON.TIME1' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.MUCA' + include 'COMMON.HAIRPIN' + j=1 + chain_beg=1 +C do i=1,nres +C write(*,*) 'initial', i,j,c(j,i) +C enddo + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then + chain_end=i + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxxsize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,chain_end + c(j,k)=c(j,k)-ireturnval*boxxsize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize + enddo + endif + chain_beg=i+1 + allareout=1 + else + if (int(c(j,i)/boxxsize).eq.0) allareout=0 + endif + enddo + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxxsize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,nres + c(j,k)=c(j,k)-ireturnval*boxxsize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize + enddo + endif +C NO JUMP +C do i=1,nres +C write(*,*) 'befor no jump', i,j,c(j,i) +C enddo + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) +C print *,'diff', difference + if (difference.gt.boxxsize/2.0) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxxsize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize + enddo + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) + if (difference.gt.boxxsize/2.0) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxxsize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxxsize + enddo + + do i=1,nres + write(*,*) 'after no jump', i,j,c(j,i) + enddo + +C NOW Y dimension + j=2 + chain_beg=1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then + chain_end=i + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxysize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,chain_end + c(j,k)=c(j,k)-ireturnval*boxysize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize + enddo + endif + chain_beg=i+1 + allareout=1 + else + if (int(c(j,i)/boxysize).eq.0) allareout=0 + endif + enddo + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxysize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,nres + c(j,k)=c(j,k)-ireturnval*boxysize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize + enddo + endif + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) + if (difference.gt.boxysize/2.0) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxysize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize + enddo + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) + if (difference.gt.boxysize/2.0) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxysize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxysize + enddo + + j=3 + chain_beg=1 + do i=1,nres-1 + if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then + chain_end=i + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxysize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,chain_end + c(j,k)=c(j,k)-ireturnval*boxzsize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize + enddo + endif + chain_beg=i+1 + allareout=1 + else + if (int(c(j,i)/boxzsize).eq.0) allareout=0 + endif + enddo + if (allareout.eq.1) then + ireturnval=int(c(j,i)/boxzsize) + if (c(j,i).le.0) ireturnval=ireturnval-1 + do k=chain_beg,nres + c(j,k)=c(j,k)-ireturnval*boxzsize + c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize + enddo + endif + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) + if (difference.gt.(boxzsize/2.0)) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxzsize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize + enddo + nojumpval=0 + do i=2,nres + if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then + difference=abs(c(j,i-1)-c(j,i)) + if (difference.gt.boxzsize/2.0) then + if (c(j,i-1).gt.c(j,i)) then + nojumpval=1 + else + nojumpval=-1 + endif + else + nojumpval=0 + endif + endif + c(j,k)=c(j,k)+nojumpval*boxzsize + c(j,k+nres)=c(j,k+nres)+nojumpval*boxzsize + enddo + + return + end diff --git a/source/wham/src-M/DIMENSIONS b/source/wham/src-M/DIMENSIONS index 154f3df..b1d6ee8 100644 --- a/source/wham/src-M/DIMENSIONS +++ b/source/wham/src-M/DIMENSIONS @@ -14,13 +14,13 @@ c parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres c parameter (maxres=250) - parameter (maxres=100) + parameter (maxres=500) C Appr. max. number of interaction sites integer maxres2 parameter (maxres2=2*maxres) C Max number of symetries integer maxsym,maxperm - parameter (maxsym=5,maxperm=120) + parameter (maxsym=12,maxperm=120) C Max. number of variables integer maxvar parameter (maxvar=4*maxres) diff --git a/source/wham/src-M/DIMENSIONS.FREE b/source/wham/src-M/DIMENSIONS.FREE index 691d9b2..e9ebe0e 100644 --- a/source/wham/src-M/DIMENSIONS.FREE +++ b/source/wham/src-M/DIMENSIONS.FREE @@ -3,8 +3,8 @@ integer MaxR,MaxT_h integer MaxSlice parameter (Max_Parm=1) - parameter (MaxQ=1,MaxQ1=MaxQ+2) - parameter(MaxR=1,MaxT_h=32) + parameter (MaxQ=4,MaxQ1=MaxQ+2) + parameter(MaxR=1,MaxT_h=36) parameter(MaxSlice=40) integer MaxN parameter (MaxN=100) diff --git a/source/wham/src-M/cinfo.f b/source/wham/src-M/cinfo.f index d497941..858334a 100644 --- a/source/wham/src-M/cinfo.f +++ b/source/wham/src-M/cinfo.f @@ -1,10 +1,10 @@ C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C -C 0 0 702 +C 0 0 703 subroutine cinfo include 'COMMON.IOUNITS' write(iout,*)'++++ Compile info ++++' - write(iout,*)'Version 0.0 build 702' - write(iout,*)'compiled Mon Dec 3 05:37:30 2012' + write(iout,*)'Version 0.0 build 703' + write(iout,*)'compiled Thu May 15 10:35:39 2014' write(iout,*)'compiled by aks255@matrix.chem.cornell.edu' write(iout,*)'OS name: Linux ' write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 ' diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F index 8e4fa74..c29d517 100644 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@ -60,6 +60,8 @@ enddo enddo enddo + write (iout,*) "indstart(me1),indend(me1)" + &,indstart(me1),indend(me1) do i=indstart(me1),indend(me1) #else do iparm=1,nParmSet @@ -71,6 +73,7 @@ enddo do i=1,ntot #endif + read(ientout,rec=i,err=101) & ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k+nres),l=1,3),k=nnt,nct), @@ -154,7 +157,10 @@ c & " kfac",kfac,"quot",quot," fT",fT & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, & wtor_d,wsccor,wbond #endif +C write (iout,*) "tuz przed energia" call etotal(energia(0),fT) +C write (iout,*) "tuz za energia" + #ifdef DEBUG write (iout,*) "Conformation",i write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) @@ -220,6 +226,7 @@ c call intout endif endif endif +C write (iout,*) "Czy tu dochodze" potE(iii+1,iparm)=energia(0) do k=1,21 enetb(k,iii+1,iparm)=energia(k) diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index a4f5fb4..810275a 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -44,11 +44,13 @@ C Gay-Berne potential (shifted LJ, angular dependence). goto 106 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence). 105 call egbv(evdw,evdw_t) +C write(iout,*) 'po elektostatyce' C C Calculate electrostatic (H-bonding) energy of the main chain. C - 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) -C + 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) +C write(iout,*) 'po eelec' + C Calculate excluded-volume interaction energy between peptide groups C and side chains. C @@ -56,8 +58,9 @@ C c c Calculate the bond-stretching energy c + call ebond(estr) -c write (iout,*) "estr",estr +C write (iout,*) "estr",estr C C Calculate the disulfide-bridge and other energy and the contributions C from other distance constraints. @@ -68,12 +71,12 @@ C C Calculate the virtual-bond-angle energy. C call ebend(ebe) -cd print *,'Bend energy finished.' +C print *,'Bend energy finished.' C C Calculate the SC local energy. C call esc(escloc) -cd print *,'SCLOC energy finished.' +C print *,'SCLOC energy finished.' C C Calculate the virtual-bond torsional energy. C @@ -1898,13 +1901,14 @@ cd enddo do i=1,nres num_cont_hb(i)=0 enddo -cd print '(a)','Enter EELEC' -cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e +C print '(a)','Enter EELEC' +C write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e do i=1,nres gel_loc_loc(i)=0.0d0 gcorr_loc(i)=0.0d0 enddo do i=iatel_s,iatel_e + if (i.eq.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & .or. itype(i+2).eq.ntyp1 & .or. itype(i-1).eq.ntyp1 @@ -1926,8 +1930,9 @@ cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=0 -c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) +C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) + if (j.eq.1) cycle if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1 & .or.itype(j+2).eq.ntyp1 & .or.itype(j-1).eq.ntyp1 @@ -2027,7 +2032,7 @@ c write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') c &'evdw1',i,j,evdwij c &,iteli,itelj,aaa,evdw1 -c write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij + write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij c write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') c & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i, c & 1.0D0/dsqrt(rrmij),evdwij,eesij, @@ -3272,7 +3277,8 @@ c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) enddo endif - + write (iout,'(a7,i5,4f7.3)') + & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff enddo estr=0.5d0*AKP*estr+estr1 c @@ -3357,6 +3363,7 @@ c write (*,'(a,i2)') 'EBEND ICG=',icg c write (iout,*) ithet_start,ithet_end do i=ithet_start,ithet_end C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) cycle if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle C Zero the energy function and its derivative at 0 or pi. @@ -3374,6 +3381,10 @@ C Zero the energy function and its derivative at 0 or pi. ichir21=isign(1,itype(i)) ichir22=isign(1,itype(i)) endif + if (i.eq.3) then + y(1)=0.0D0 + y(2)=0.0D0 + else if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF @@ -3390,6 +3401,7 @@ C Zero the energy function and its derivative at 0 or pi. y(1)=0.0D0 y(2)=0.0D0 endif + endif if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) @@ -3457,6 +3469,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai +c write (iout,'(a6,i5,0pf7.3,f7.3,i5)') +c & 'ebend',i,ethetai,theta(i),itype(i) c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i), c & rad2deg*phii,rad2deg*phii1,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 @@ -3605,7 +3619,9 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end +C if (i.eq.2) cycle C if (itype(i-1).eq.ntyp1) cycle + if (i.le.2) cycle if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1 & .or.itype(i).eq.ntyp1) cycle if (iabs(itype(i+1)).eq.20) iblock=2 @@ -3619,6 +3635,14 @@ C if (itype(i-1).eq.ntyp1) cycle coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo + if (i.eq.3) then + phii=0.0d0 + ityp1=nthetyp+1 + do k=1,nsingle + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + else if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) @@ -3639,6 +3663,7 @@ C if (itype(i-1).eq.ntyp1) cycle sinph1(k)=0.0d0 enddo endif + endif if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) @@ -3799,14 +3824,14 @@ C ALPHA and OMEGA. common /sccalc/ time11,time12,time112,theti,it,nlobit delta=0.02d0*pi escloc=0.0D0 -c write (iout,'(a)') 'ESC' +C write (iout,*) 'ESC' do i=loc_start,loc_end it=itype(i) if (it.eq.ntyp1) cycle if (it.eq.10) goto 1 nlobit=nlob(iabs(it)) c print *,'i=',i,' it=',it,' nlobit=',nlobit -c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad +C write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad theti=theta(i+1)-pipol x(1)=dtan(theti) x(2)=alph(i) @@ -3842,8 +3867,8 @@ c write (iout,*) "i",i," x",x(1),x(2),x(3) dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) -c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd + write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, + & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi c escloci=esclocbi c write (iout,*) escloci @@ -3877,15 +3902,17 @@ c write (iout,*) escloci enddo dersc(2)=dersc(2)+ssd*(escloci-esclocbi) c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, -c & esclocbi,ss,ssd +c & esclocbi,ss,ssd escloci=ss*escloci+(1.0d0-ss)*esclocbi -c write (iout,*) escloci +C write (iout,*) 'i=',i, escloci else call enesc(x,escloci,dersc,ddummy,.false.) endif escloc=escloc+escloci -c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc +C write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc + write (iout,'(a6,i5,0pf7.3)') + & 'escloc',i,escloci gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & wscloc*dersc(1) @@ -4575,6 +4602,7 @@ C Set lprn=.true. for debugging c lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end + if (i.le.2) cycle if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle C if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 @@ -4678,6 +4706,7 @@ C Set lprn=.true. for debugging c lprn=.true. etors_d=0.0D0 do i=iphi_start,iphi_end-1 + if (i.le.3) cycle C if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 C & .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or. diff --git a/source/wham/src-M/include_unres/COMMON.INTERACT b/source/wham/src-M/include_unres/COMMON.INTERACT index 650af0d..d38c470 100644 --- a/source/wham/src-M/include_unres/COMMON.INTERACT +++ b/source/wham/src-M/include_unres/COMMON.INTERACT @@ -20,7 +20,10 @@ C 12/1/95 Array EPS included in the COMMON block. & eps_scp(ntyp,2),rscp(ntyp,2),eps_orig(ntyp,ntyp) c 12/5/03 modified 09/18/03 Bond stretching parameters. double precision vbldp0,vbldsc0,akp,aksc,abond0,distchainmax + &,vbldpDUM integer nbondterm common /stretch/ vbldp0,vbldsc0(maxbondterm,ntyp),akp, & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) + &,vbldpDUM + diff --git a/source/wham/src-M/initialize_p.F b/source/wham/src-M/initialize_p.F index 9d159a9..d1cc94b 100644 --- a/source/wham/src-M/initialize_p.F +++ b/source/wham/src-M/initialize_p.F @@ -537,7 +537,7 @@ C Partition local interactions #else loc_start=2 loc_end=nres-1 - ithet_start=3 + ithet_start=3 ithet_end=nres iphi_start=nnt+3 iphi_end=nct