From: Adam Sieradzan Date: Wed, 11 Feb 2015 11:23:12 +0000 (+0100) Subject: intro of SCCOR into src_MIN X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=5c15a6184443adadbb087d7badecb43cd04a77db intro of SCCOR into src_MIN --- diff --git a/source/unres/src_CSA/CMakeLists.txt b/source/unres/src_CSA/CMakeLists.txt index fdea19f..b297d98 100644 --- a/source/unres/src_CSA/CMakeLists.txt +++ b/source/unres/src_CSA/CMakeLists.txt @@ -146,7 +146,7 @@ elseif(UNRES_CSA_FF STREQUAL "ALPHABETA") elseif(UNRES_CSA_FF STREQUAL "CASP5") set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DMOMENT -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" ) elseif(UNRES_CSA_FF STREQUAL "3P") - set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" ) + set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" ) elseif(UNRES_CSA_FF STREQUAL "4P") set(CPPFLAGS "UNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB" ) endif(UNRES_CSA_FF STREQUAL "CASP3") diff --git a/source/unres/src_MD/CMakeLists.txt b/source/unres/src_MD/CMakeLists.txt index ea7776c..fadf9ca 100644 --- a/source/unres/src_MD/CMakeLists.txt +++ b/source/unres/src_MD/CMakeLists.txt @@ -166,14 +166,14 @@ set_property(SOURCE ${UNRES_MD_SRC3} PROPERTY COMPILE_FLAGS ${FFLAGS3} ) #========================================= if(UNRES_MD_FF STREQUAL "GAB" ) # set preprocesor flags - set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC" ) + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB") #========================================= # Settings for E0LL2Y force field #========================================= elseif(UNRES_MD_FF STREQUAL "E0LL2Y") # set preprocesor flags - set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0" ) + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB" ) endif(UNRES_MD_FF STREQUAL "GAB") #========================================= diff --git a/source/unres/src_MD/readpdb.F b/source/unres/src_MD/readpdb.F index 48e0abd..80c3fa3 100644 --- a/source/unres/src_MD/readpdb.F +++ b/source/unres/src_MD/readpdb.F @@ -252,7 +252,20 @@ ctest stop vbld(i+1)=dist(i,i+1) vbld_inv(i+1)=1.0d0/vbld(i+1) if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1) - if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + if (i.gt.2) then + if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1) + if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then + tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres) + endif + if (itype(i-1).ne.10) then + tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1) + omicron(1,i)=alpha(i-2,i-1,i-1+nres) + omicron(2,i)=alpha(i-1+nres,i-1,i) + endif + if (itype(i).ne.10) then + tauangle(2,i+1)=beta(i-2,i-1,i,i+nres) + endif + endif enddo c if (unres_pdb) then c if (itype(1).eq.21) then diff --git a/source/unres/src_MIN/COMMON.CHAIN b/source/unres/src_MIN/COMMON.CHAIN index f7a8a1d..6e19f8d 100644 --- a/source/unres/src_MIN/COMMON.CHAIN +++ b/source/unres/src_MIN/COMMON.CHAIN @@ -1,9 +1,10 @@ integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc, & nres0,nstart_seq double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, - & prod,rt,dc_work,cref,crefjlee + & prod,rt,dc_work,cref,crefjlee,dc_norm2 common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), + & dc_norm2(3,0:maxres2), & dc_work(MAXRES6),nres,nres0 common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres), & rt(3,3,maxres) diff --git a/source/unres/src_MIN/COMMON.LOCAL b/source/unres/src_MIN/COMMON.LOCAL index 837a7a3..2c9ea5f 100644 --- a/source/unres/src_MIN/COMMON.LOCAL +++ b/source/unres/src_MIN/COMMON.LOCAL @@ -32,7 +32,7 @@ C Virtual-bond lenghts & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end, & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end, & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start, - & iint_end,iphi1_start,iphi1_end, + & iint_end,iphi1_start,iphi1_end,itau_start,itau_end, & ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1), & ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1), & iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1), @@ -46,7 +46,7 @@ C Virtual-bond lenghts & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end, & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start, & iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ, - & ivec_count,iset_displ, + & ivec_count,iset_displ,itau_start,itau_end, & iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count, & iphi_displ,iphi_count,iphi1_displ,iphi1_count C Inverses of the actual virtual bond lengths diff --git a/source/unres/src_MIN/COMMON.SCCOR b/source/unres/src_MIN/COMMON.SCCOR index a28f621..e29cb4c 100644 --- a/source/unres/src_MIN/COMMON.SCCOR +++ b/source/unres/src_MIN/COMMON.SCCOR @@ -1,6 +1,18 @@ -C Parameters of the SCCOR term - double precision v1sccor,v2sccor - integer nterm_sccor - common/sccor/v1sccor(maxterm_sccor,20,20), - & v2sccor(maxterm_sccor,20,20), - & nterm_sccor +cc Parameters of the SCCOR term + double precision v1sccor,v2sccor,vlor1sccor, + & vlor2sccor,vlor3sccor,gloc_sc, + & dcostau,dsintau,dtauangle,dcosomicron, + & domicron + integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor + common/sccor/v1sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp), + & v2sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp), + & v0sccor(maxterm_sccor,-ntyp:ntyp), + & nterm_sccor(-ntyp:ntyp,-ntyp:ntyp),isccortyp(-ntyp:ntyp), + & nsccortyp, + & nlor_sccor(-ntyp:ntyp,-ntyp:ntyp), + & vlor1sccor(maxterm_sccor,20,20), + & vlor2sccor(maxterm_sccor,20,20), + & vlor3sccor(maxterm_sccor,20,20),gloc_sc(3,0:maxres2,10), + & dcostau(3,3,3,maxres2),dsintau(3,3,3,maxres2), + & dtauangle(3,3,3,maxres2),dcosomicron(3,3,3,maxres2), + & domicron(3,3,3,maxres2) diff --git a/source/unres/src_MIN/COMMON.VAR b/source/unres/src_MIN/COMMON.VAR index 71158b8..fd0626a 100644 --- a/source/unres/src_MIN/COMMON.VAR +++ b/source/unres/src_MIN/COMMON.VAR @@ -3,11 +3,13 @@ C Store the geometric variables in the following COMMON block. & mask_theta,mask_phi,mask_side double precision theta,phi,alph,omeg,varsave,esave,varall,vbld, & thetaref,phiref,costtab,sinttab,cost2tab,sint2tab, + & tauangle,omicron, & xxtab,yytab,zztab,xxref,yyref,zzref common /var/ theta(maxres),phi(maxres),alph(maxres),omeg(maxres), & vbld(2*maxres),thetaref(maxres),phiref(maxres), & costtab(maxres), sinttab(maxres), cost2tab(maxres), & sint2tab(maxres),xxtab(maxres),yytab(maxres), + & omicron(2,maxres),tauangle(3,maxres), & zztab(maxres),xxref(maxres),yyref(maxres),zzref(maxres), & ialph(maxres,2),ivar(4*maxres2),ntheta,nphi,nside,nvar C Store the angles and variables corresponding to old conformations (for use diff --git a/source/unres/src_MIN/checkder_p.F b/source/unres/src_MIN/checkder_p.F index 67465bb..a50bc2a 100644 --- a/source/unres/src_MIN/checkder_p.F +++ b/source/unres/src_MIN/checkder_p.F @@ -514,7 +514,20 @@ c------------------------------------------------------------------------- & +(c(j,i+1)-c(j,i))/dnorm2) enddo be=0.0D0 - if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + if (i.gt.2) then + if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1) + if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then + tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres) + endif + if (itype(i-1).ne.10) then + tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1) + omicron(1,i)=alpha(i-2,i-1,i-1+nres) + omicron(2,i)=alpha(i-1+nres,i-1,i) + endif + if (itype(i).ne.10) then + tauangle(2,i+1)=beta(i-2,i-1,i,i+nres) + endif + endif omeg(i)=beta(nres+i,i,maxres2,i+1) alph(i)=alpha(nres+i,i,maxres2) theta(i+1)=alpha(i-1,i,i+1) diff --git a/source/unres/src_MIN/energy_p_new_barrier.F b/source/unres/src_MIN/energy_p_new_barrier.F index 96df440..3f19648 100644 --- a/source/unres/src_MIN/energy_p_new_barrier.F +++ b/source/unres/src_MIN/energy_p_new_barrier.F @@ -5780,29 +5780,53 @@ c amino-acid residues. include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging +C Set lprn=.true. for debugging lprn=.false. c lprn=.true. -c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor +c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 - do i=iphi_start,iphi_end + do i=itau_start,itau_end esccor_ii=0.0D0 - itori=itype(i-2) - itori1=itype(i-1) + isccori=isccortyp(itype(i-2)) + isccori1=isccortyp(itype(i-1)) +c write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1) phii=phi(i) + do intertyp=1,3 !intertyp +cc Added 09 May 2012 (Adasko) +cc Intertyp means interaction type of backbone mainchain correlation: +c 1 = SC...Ca...Ca...Ca +c 2 = Ca...Ca...Ca...SC +c 3 = SC...Ca...Ca...SCi gloci=0.0D0 - do j=1,nterm_sccor - v1ij=v1sccor(j,itori,itori1) - v2ij=v2sccor(j,itori,itori1) - cosphi=dcos(j*phii) - sinphi=dsin(j*phii) + if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. + & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-1).eq.ntyp1))) + & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) + & .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) + & .or.(itype(i).eq.ntyp1))) + & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. + & (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + & (itype(i-3).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) + & cycle + do j=1,nterm_sccor(isccori,isccori1) + v1ij=v1sccor(j,intertyp,isccori,isccori1) + v2ij=v2sccor(j,intertyp,isccori,isccori1) + cosphi=dcos(j*tauangle(intertyp,i)) + sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo +c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp + gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') - & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, - & (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6) - gsccor_loc(i-3)=gsccor_loc(i-3)+gloci + & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1, + & (v1sccor(j,intertyp,isccori,isccori1),j=1,6) + & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6) +C gsccor_loc(i-3)=gsccor_loc(i-3)+gloci + enddo !intertyp enddo return end diff --git a/source/unres/src_MIN/gradient_p.F b/source/unres/src_MIN/gradient_p.F index 25d1b12..569120f 100644 --- a/source/unres/src_MIN/gradient_p.F +++ b/source/unres/src_MIN/gradient_p.F @@ -337,6 +337,7 @@ C------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.VAR' include 'COMMON.MD_' + include 'COMMON.SCCOR' C C Initialize Cartesian-coordinate gradient C @@ -374,6 +375,9 @@ C gradx(j,i,icg)=0.0d0 gscloc(j,i)=0.0d0 gsclocx(j,i)=0.0d0 + do intertyp=1,3 + gloc_sc(intertyp,i,icg)=0.0d0 + enddo enddo enddo C diff --git a/source/unres/src_MIN/initialize_p.F b/source/unres/src_MIN/initialize_p.F index 8cc2278..8cca656 100644 --- a/source/unres/src_MIN/initialize_p.F +++ b/source/unres/src_MIN/initialize_p.F @@ -564,6 +564,9 @@ C Partition local interactions iphi_end=iturn3_end+2 iturn3_start=iturn3_start-1 iturn3_end=iturn3_end-1 + call int_bounds(nres-3,itau_start,itau_end) + itau_start=itau_start+3 + itau_end=itau_end+3 call int_bounds(nres-3,iphi1_start,iphi1_end) iphi1_start=iphi1_start+3 iphi1_end=iphi1_end+3 @@ -1090,6 +1093,8 @@ c write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 idihconstr_end=ndih_constr iphid_start=iphi_start iphid_end=iphi_end-1 + itau_start=4 + itau_end=nres ibond_start=2 ibond_end=nres-1 ibondp_start=nnt+1 diff --git a/source/unres/src_MIN/int_to_cart.f b/source/unres/src_MIN/int_to_cart.f index 97324ec..db8aeb0 100644 --- a/source/unres/src_MIN/int_to_cart.f +++ b/source/unres/src_MIN/int_to_cart.f @@ -13,9 +13,9 @@ c------------------------------------------------------------- include 'COMMON.INTERACT' include 'COMMON.MD_' include 'COMMON.IOUNITS' - + include 'COMMON.SCCOR' c calculating dE/ddc1 - if (nres.lt.3) return + if (nres.lt.3) go to 18 do j=1,3 gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) & +gloc(nres-2,icg)*dtheta(j,1,3) @@ -112,7 +112,162 @@ c The side-chain vector derivatives & +gloc(ialph(i,1)+nside,icg)*domega(j,3,i) enddo endif - enddo + enddo +c---------------------------------------------------------------------- +C INTERTYP=1 SC...Ca...Ca...Ca +C INTERTYP=2 Ca...Ca...Ca...SC +C INTERTYP=3 SC...Ca...Ca...SC +c calculating dE/ddc1 + 18 continue +c do i=1,nres +c gloc(i,icg)=0.0D0 +c write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg) +c enddo + if (nres.lt.2) return + if ((nres.lt.3).and.(itype(1).eq.10)) return + if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then + do j=1,3 +cc Derviative was calculated for oposite vector of side chain therefore +c there is "-" sign before gloc_sc + gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)* + & dtauangle(j,1,1,3) + gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* + & dtauangle(j,1,2,3) + if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then + gxcart(j,1)= gxcart(j,1) + & -gloc_sc(3,0,icg)*dtauangle(j,3,1,3) + gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* + & dtauangle(j,3,2,3) + endif + enddo + endif + if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) + & then + do j=1,3 + gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4) + enddo + endif +c As potetnial DO NOT depend on omicron anlge their derivative is +c ommited +c & +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3) + +c Calculating the remainder of dE/ddc2 + do j=1,3 + if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then + if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ + & gloc_sc(3,0,icg)*dtauangle(j,3,3,3) + if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) + & then + gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4) +cc the - above is due to different vector direction + gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4) + endif + if (nres.gt.3) then + gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4) +cc the - above is due to different vector direction + gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4) +c write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart" +c write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx" + endif + endif + if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then + gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3) +c write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3) + endif + if ((itype(3).ne.10).and.(nres.ge.3)) then + gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4) +c write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4) + endif + if ((itype(4).ne.10).and.(nres.ge.4)) then + gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5) +c write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5) + endif + +c write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2" + enddo +c If there are more than five residues + if(nres.ge.5) then + do i=3,nres-2 + do j=1,3 +c write(iout,*) "before", gcart(j,i) + if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) + & *dtauangle(j,2,3,i+1) + & -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2) + gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg) + & *dtauangle(j,1,2,i+2) +c write(iout,*) "new",j,i, +c & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2) + if (itype(i-1).ne.10) then + gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) + &*dtauangle(j,3,3,i+1) + endif + if (itype(i+1).ne.10) then + gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) + &*dtauangle(j,3,1,i+2) + gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) + &*dtauangle(j,3,2,i+2) + endif + endif + if (itype(i-1).ne.10) then + gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* + & dtauangle(j,1,3,i+1) + endif + if (itype(i+1).ne.10) then + gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* + & dtauangle(j,2,2,i+2) +c write(iout,*) "numer",i,gloc_sc(2,i-1,icg), +c & dtauangle(j,2,2,i+2) + endif + if (itype(i+2).ne.10) then + gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* + & dtauangle(j,2,1,i+3) + endif + enddo + enddo + endif +c Setting dE/ddnres-1 + if(nres.ge.4) then + do j=1,3 + if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then + gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) + & *dtauangle(j,2,3,nres) +c write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg), +c & dtauangle(j,2,3,nres), gxcart(j,nres-1) + if (itype(nres-2).ne.10) then + gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) + & *dtauangle(j,3,3,nres) + endif + if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then + gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) + & *dtauangle(j,3,1,nres+1) + gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) + & *dtauangle(j,3,2,nres+1) + endif + endif + if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then + gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* + & dtauangle(j,1,3,nres) + endif + if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then + gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* + & dtauangle(j,2,2,nres+1) +c write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg), +c & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres) + endif + enddo + endif +c Settind dE/ddnres + if ((nres.ge.3).and.(itype(nres).ne.10).and. + & (itype(nres).ne.ntyp1))then + do j=1,3 + gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) + & *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) + & *dtauangle(j,2,3,nres+1) + enddo + endif +c The side-chain vector derivatives + return end diff --git a/source/unres/src_MIN/intcartderiv.F b/source/unres/src_MIN/intcartderiv.F index 5fea875..2c865e6 100644 --- a/source/unres/src_MIN/intcartderiv.F +++ b/source/unres/src_MIN/intcartderiv.F @@ -12,6 +12,7 @@ include 'COMMON.DERIV' include 'COMMON.IOUNITS' include 'COMMON.LOCAL' + include 'COMMON.SCCOR' double precision dcostheta(3,2,maxres), & dcosphi(3,3,maxres),dsinphi(3,3,maxres), & dcosalpha(3,3,maxres),dcosomega(3,3,maxres), @@ -44,7 +45,41 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i) dtheta(j,2,i)=-1/sint*dcostheta(j,2,i) enddo enddo - +#if defined(MPI) && defined(PARINTDER) +c We need dtheta(:,:,i-1) to compute dphi(:,:,i) + do i=max0(ithet_start-1,3),ithet_end +#else + do i=3,nres +#endif + if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then + cost1=dcos(omicron(1,i)) + sint1=sqrt(1-cost1*cost1) + cost2=dcos(omicron(2,i)) + sint2=sqrt(1-cost2*cost2) + do j=1,3 +CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1) + dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ + & cost1*dc_norm(j,i-2))/ + & vbld(i-1) + domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i) + dcosomicron(j,1,2,i)=-(dc_norm(j,i-2) + & +cost1*(dc_norm(j,i-1+nres)))/ + & vbld(i-1+nres) + domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i) +CC Calculate derivative over second omicron Sci-1,Cai-1 Cai +CC Looks messy but better than if in loop + dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres) + & +cost2*dc_norm(j,i-1))/ + & vbld(i) + domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i) + dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) + & +cost2*(-dc_norm(j,i-1+nres)))/ + & vbld(i-1+nres) +c write(iout,*) "vbld", i,itype(i),vbld(i-1+nres) + domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i) + enddo + endif + enddo c Derivatives of phi: c If phi is 0 or 180 degrees, then the formulas c have to be derived by power series expansion of the @@ -109,6 +144,227 @@ c Obtaining the gamma derivatives from cosine derivative enddo endif enddo +Calculate derivative of Tauangle +#ifdef PARINTDER + do i=itau_start,itau_end +#else + do i=3,nres +#endif + if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle +c if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or. +c & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle +cc dtauangle(j,intertyp,dervityp,residue number) +cc INTERTYP=1 SC...Ca...Ca..Ca +c the conventional case + sint=dsin(theta(i)) + sint1=dsin(omicron(2,i-1)) + sing=dsin(tauangle(1,i)) + cost=dcos(theta(i)) + cost1=dcos(omicron(2,i-1)) + cosg=dcos(tauangle(1,i)) + do j=1,3 + dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres) +cc write(iout,*) dc_norm2(j,i-2+nres),"dcnorm" + enddo + scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1)) + fac0=1.0d0/(sint1*sint) + fac1=cost*fac0 + fac2=cost1*fac0 + fac3=cosg*cost1/(sint1*sint1) + fac4=cosg*cost/(sint*sint) +cc write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4 +c Obtaining the gamma derivatives from sine derivative + if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. + & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. + & tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then + call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1) + call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2) + call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3) + do j=1,3 + ctgt=cost/sint + ctgt1=cost1/sint1 + cosg_inv=1.0d0/cosg + dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1) + &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres))) + & *vbld_inv(i-2+nres) + dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i) + dsintau(j,1,2,i)= + & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i)) + & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1) +c write(iout,*) "dsintau", dsintau(j,1,2,i) + dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i) +c Bug fixed 3/24/05 (AL) + dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i) + & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i) +c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) + dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i) + enddo +c Obtaining the gamma derivatives from cosine derivative + else + do j=1,3 + dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3* + & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp* + & (dc_norm2(j,i-2+nres)))/vbld(i-2+nres) + dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i) + dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2* + & dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4* + & dcostheta(j,1,i) + dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i) + dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4* + & dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp* + & dc_norm(j,i-1))/vbld(i) + dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i) +c write (iout,*) "else",i + enddo + endif +c do k=1,3 +c write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3) +c enddo + enddo +CC Second case Ca...Ca...Ca...SC +#ifdef PARINTDER + do i=itau_start,itau_end +#else + do i=4,nres +#endif + if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. + & (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle +c the conventional case + sint=dsin(omicron(1,i)) + sint1=dsin(theta(i-1)) + sing=dsin(tauangle(2,i)) + cost=dcos(omicron(1,i)) + cost1=dcos(theta(i-1)) + cosg=dcos(tauangle(2,i)) +c do j=1,3 +c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) +c enddo + scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres)) + fac0=1.0d0/(sint1*sint) + fac1=cost*fac0 + fac2=cost1*fac0 + fac3=cosg*cost1/(sint1*sint1) + fac4=cosg*cost/(sint*sint) +c Obtaining the gamma derivatives from sine derivative + if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or. + & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or. + & tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then + call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1) + call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2) + call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3) + do j=1,3 + ctgt=cost/sint + ctgt1=cost1/sint1 + cosg_inv=1.0d0/cosg + dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1) + & +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2) +c write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1), +c &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)" + dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i) + dsintau(j,2,2,i)= + & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i)) + & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1) +c write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1), +c & sing*ctgt*domicron(j,1,2,i), +c & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1) + dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i) +c Bug fixed 3/24/05 (AL) + dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i) + & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres) +c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) + dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i) + enddo +c Obtaining the gamma derivatives from cosine derivative + else + do j=1,3 + dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3* + & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp* + & dc_norm(j,i-3))/vbld(i-2) + dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i) + dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2* + & dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4* + & dcosomicron(j,1,1,i) + dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i) + dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4* + & dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp* + & dc_norm(j,i-1+nres))/vbld(i-1+nres) + dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i) +c write(iout,*) i,j,"else", dtauangle(j,2,3,i) + enddo + endif + enddo + +CCC third case SC...Ca...Ca...SC +#ifdef PARINTDER + + do i=itau_start,itau_end +#else + do i=3,nres +#endif +c the conventional case + if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. + &(itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle + sint=dsin(omicron(1,i)) + sint1=dsin(omicron(2,i-1)) + sing=dsin(tauangle(3,i)) + cost=dcos(omicron(1,i)) + cost1=dcos(omicron(2,i-1)) + cosg=dcos(tauangle(3,i)) + do j=1,3 + dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres) +c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres) + enddo + scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres)) + fac0=1.0d0/(sint1*sint) + fac1=cost*fac0 + fac2=cost1*fac0 + fac3=cosg*cost1/(sint1*sint1) + fac4=cosg*cost/(sint*sint) +c Obtaining the gamma derivatives from sine derivative + if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. + & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. + & tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then + call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1) + call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2) + call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3) + do j=1,3 + ctgt=cost/sint + ctgt1=cost1/sint1 + cosg_inv=1.0d0/cosg + dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1) + & -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres)) + & *vbld_inv(i-2+nres) + dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i) + dsintau(j,3,2,i)= + & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i)) + & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1) + dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i) +c Bug fixed 3/24/05 (AL) + dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i) + & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres)) + & *vbld_inv(i-1+nres) +c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1) + dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i) + enddo +c Obtaining the gamma derivatives from cosine derivative + else + do j=1,3 + dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3* + & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp* + & dc_norm2(j,i-2+nres))/vbld(i-2+nres) + dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i) + dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2* + & dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4* + & dcosomicron(j,1,1,i) + dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i) + dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4* + & dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp* + & dc_norm(j,i-1+nres))/vbld(i-1+nres) + dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i) +c write(iout,*) "else",i + enddo + endif + enddo #ifdef CRYST_SC c Derivatives of side-chain angles alpha and omega #if defined(MPI) && defined(PARINTDER) diff --git a/source/unres/src_MIN/parmread.F b/source/unres/src_MIN/parmread.F index f95bbb0..20a2319 100644 --- a/source/unres/src_MIN/parmread.F +++ b/source/unres/src_MIN/parmread.F @@ -557,24 +557,126 @@ C C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local C correlation energies. C - read (isccor,*,end=119,err=119) nterm_sccor - do i=1,20 - do j=1,20 - read (isccor,'(a)') - do k=1,nterm_sccor - read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j), - & v2sccor(k,i,j) + read (isccor,*,end=119,err=119) nsccortyp +#ifdef SCCORPDB + read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp) +C For D-aminoacid uncomment +C do i=-ntyp,-1 +C isccortyp(i)=-isccortyp(-i) +C enddo + iscprol=isccortyp(20) +c write (iout,*) 'ntortyp',ntortyp + maxinter=3 +cc maxinter is maximum interaction sites + do l=1,maxinter + do i=1,nsccortyp + do j=1,nsccortyp + read (isccor,*,end=119,err=119) + &nterm_sccor(i,j),nlor_sccor(i,j) + v0ijsccor=0.0d0 + si=-1.0d0 +C nterm_sccor(-i,j)=nterm_sccor(i,j) +C nterm_sccor(-i,-j)=nterm_sccor(i,j) +C nterm_sccor(i,-j)=nterm_sccor(i,j) + do k=1,nterm_sccor(i,j) + read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) +C if (j.eq.iscprol) then +C if (i.eq.isccortyp(10)) then +C v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) +C v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) +C else +C v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 +C & +v2sccor(k,l,i,j)*dsqrt(0.75d0) +C v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 +C & +v1sccor(k,l,i,j)*dsqrt(0.75d0) +C v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) +C v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) +C v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) +C v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) +C endif +CC else +C if (i.eq.isccortyp(10)) then +C v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) +C v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) +C else +C if (j.eq.isccortyp(10)) then +C v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) +C v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) +C else +C v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) +C v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) +C v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) +C v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) +C v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) +C v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) +C endif +C endif +C endif + v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + si=-si + enddo + do k=1,nlor_sccor(i,j) + read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j), + & vlor2sccor(k,i,j),vlor3sccor(k,i,j) + v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ + &(1+vlor3sccor(k,i,j)**2) + enddo + v0sccor(i,j)=v0ijsccor + enddo + enddo + enddo + close (isccor) +#else + read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp) +c write (iout,*) 'ntortyp',ntortyp + maxinter=3 +cc maxinter is maximum interaction sites + do l=1,maxinter + do i=1,nsccortyp + do j=1,nsccortyp + read (isccor,*,end=113,err=113) + & nterm_sccor(i,j),nlor_sccor(i,j) + v0ijsccor=0.0d0 + si=-1.0d0 + + do k=1,nterm_sccor(i,j) + read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j) + & ,v2sccor(k,l,i,j) + v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j) + si=-si enddo + do k=1,nlor_sccor(i,j) + read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j), + & vlor2sccor(k,i,j),vlor3sccor(k,i,j) + v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ + &(1+vlor3sccor(k,i,j)**2) + enddo + v0sccor(i,j)=v0ijsccor enddo enddo + enddo close (isccor) + +#endif if (lprint) then - write (iout,'(/a/)') 'Torsional constants of SCCORR:' - do i=1,20 - do j=1,20 + write (iout,'(/a/)') 'Torsional constants of SCCORR:' +C do i=1,20 +C do j=1,20 +C write (iout,'(/a/)') 'Torsional constants:' + do i=1,nsccortyp + do j=1,nsccortyp write (iout,*) 'ityp',i,' jtyp',j - do k=1,nterm_sccor - write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j) +C do k=1,nterm_sccor +C write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j) +C write (iout,*) 'Fourier constants' + do k=1,nterm_sccor(i,j) + write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j) + enddo + write (iout,*) 'Lorenz constants' + do k=1,nlor_sccor(i,j) + write (iout,'(3(1pe15.5))') + & vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j) enddo enddo enddo diff --git a/source/unres/src_MIN/readpdb.F b/source/unres/src_MIN/readpdb.F index eb4ba3f..8dec8a8 100644 --- a/source/unres/src_MIN/readpdb.F +++ b/source/unres/src_MIN/readpdb.F @@ -251,6 +251,7 @@ ctest stop vbld_inv(i+1)=1.0d0/vbld(i+1) if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1) if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + enddo c if (unres_pdb) then c if (itype(1).eq.21) then