From: Adam Kazimierz Sieradzan Date: Thu, 4 Oct 2012 14:52:15 +0000 (-0400) Subject: Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=2acc9917f7452c3a59a242e6033bea16f846cbfe Wprowadznie potencjalow SC-COR do multichain oraz ich pseudosymetrii dla D aminokwasow --- diff --git a/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe b/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe index 0cc9b64..71157fe 100755 Binary files a/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe and b/bin/unres/MD-M/unres_Tc_procor_oldparm_em64-D-symetr.exe differ diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index ea1c17a..71db8f7 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -1,9 +1,10 @@ integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc, & nres0,nstart_seq,chain_length,iprzes,tabperm,nperm double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, - & prod,rt,dc_work,cref,crefjlee,chain_rep + & prod,rt,dc_work,cref,crefjlee,chain_rep,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_MD-M/COMMON.LOCAL b/source/unres/src_MD-M/COMMON.LOCAL index e1e13f4..c8d686f 100644 --- a/source/unres/src_MD-M/COMMON.LOCAL +++ b/source/unres/src_MD-M/COMMON.LOCAL @@ -34,7 +34,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), @@ -48,7 +48,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_MD-M/COMMON.SCCOR b/source/unres/src_MD-M/COMMON.SCCOR index 5217de7..e29cb4c 100644 --- a/source/unres/src_MD-M/COMMON.SCCOR +++ b/source/unres/src_MD-M/COMMON.SCCOR @@ -1,6 +1,18 @@ -C Parameters of the SCCOR term - double precision v1sccor,v2sccor - integer nterm_sccor - common/torsion/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_MD-M/COMMON.VAR b/source/unres/src_MD-M/COMMON.VAR index 71158b8..1ab0a16 100644 --- a/source/unres/src_MD-M/COMMON.VAR +++ b/source/unres/src_MD-M/COMMON.VAR @@ -3,10 +3,12 @@ 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), + & omicron(2,maxres),tauangle(3,maxres), & sint2tab(maxres),xxtab(maxres),yytab(maxres), & zztab(maxres),xxref(maxres),yyref(maxres),zzref(maxres), & ialph(maxres,2),ivar(4*maxres2),ntheta,nphi,nside,nvar diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index d9992af..5131f4e 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -59,7 +59,7 @@ C virtual-bond angle bending potentials & mmaxtheterm=maxtheterm) c Max number of torsional terms in SCCOR integer maxterm_sccor - parameter (maxterm_sccor=3) + parameter (maxterm_sccor=6) C Max. number of lobes in SC distribution integer maxlob parameter (maxlob=4) diff --git a/source/unres/src_MD-M/MD.F b/source/unres/src_MD-M/MD.F index a28d660..f65a527 100644 --- a/source/unres/src_MD-M/MD.F +++ b/source/unres/src_MD-M/MD.F @@ -1628,7 +1628,7 @@ c if the friction coefficients do not depend on surface area stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct - stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) + stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i)))) enddo endif c Open the pdb file for snapshotshots diff --git a/source/unres/src_MD-M/Makefile b/source/unres/src_MD-M/Makefile index e09600e..35006b3 100644 --- a/source/unres/src_MD-M/Makefile +++ b/source/unres/src_MD-M/Makefile @@ -1,7 +1,7 @@ CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI \ -DPGI -DSPLITELE -DISNAN -DAMD64 \ -DPROCOR -DLANG0 \ - -DCRYST_BOND -DCRYST_THETA -DCRYST_SC + -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB ## -DPROCOR ## -DMOMENT #-DCO_BIAS @@ -14,8 +14,9 @@ INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh #FC= /usr/local/opt/intel/compiler60/ia32/bin/ifc FC= ifort -OPT = -O3 -ip -w -#OPT = -g -CB +#OPT = -O3 -ip -w +OPT = -g -CB +#OPT = -g CFLAGS = -DSGI -c FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index 26854e6..0539e48 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -513,7 +513,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_MD-M/cinfo.f b/source/unres/src_MD-M/cinfo.f index f46c7f9..4e7d75d 100644 --- a/source/unres/src_MD-M/cinfo.f +++ b/source/unres/src_MD-M/cinfo.f @@ -1,10 +1,10 @@ C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C -C 2 3 3404 +C 2 3 3438 subroutine cinfo include 'COMMON.IOUNITS' write(iout,*)'++++ Compile info ++++' - write(iout,*)'Version 2.3 build 3404' - write(iout,*)'compiled Thu Aug 16 11:51:48 2012' + write(iout,*)'Version 2.3 build 3438' + write(iout,*)'compiled Thu Oct 4 09:08:03 2012' 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 ' @@ -14,10 +14,10 @@ C 2 3 3404 write(iout,*)'CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI \\' write(iout,*)' -DPGI -DSPLITELE -DISNAN -DAMD64 \\' write(iout,*)' -DPROCOR -DLANG0 \\' - write(iout,*)' -DCRYST_BOND -DCRYST_THETA -DCRYST_SC ' + write(iout,*)' -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORP...' write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...' write(iout,*)'FC= ifort' - write(iout,*)'OPT = -O3 -ip -w' + write(iout,*)'OPT = -g -CB' write(iout,*)'CFLAGS = -DSGI -c' write(iout,*)'FFLAGS = -c ${OPT} -I$(INSTALL_DIR)/include' write(iout,*)'FFLAGS1 = -c -w -g -d2 -CA -CB -I$(INSTALL_DIR)...' diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 0ce1a7b..38c5306 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -434,7 +434,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include 'mpif.h' double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres) + & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres) #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' @@ -447,6 +447,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.TIME1' include 'COMMON.MAXGRAD' + include 'COMMON.SCCOR' #ifdef TIMING time01=MPI_Wtime() #endif @@ -689,7 +690,6 @@ c enddo & +wturn3*gel_loc_turn3(i) & +wturn6*gel_loc_turn6(i) & +wel_loc*gel_loc_loc(i) - & +wsccor*gsccor_loc(i) enddo #ifdef DEBUG write (iout,*) "gloc after adding corr" @@ -708,6 +708,21 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo +#define DEBUG +#ifdef DEBUG + write (iout,*) "gloc_sc before reduce" + do i=1,nres + do j=1,1 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif +#undef DEBUG + do i=1,nres + do j=1,3 + gloc_scbuf(j,i)=gloc_sc(j,i,icg) + enddo + enddo time00=MPI_Wtime() call MPI_Barrier(FG_COMM,IERR) time_barrier_g=time_barrier_g+MPI_Wtime()-time00 @@ -719,6 +734,19 @@ c enddo call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) time_reduce=time_reduce+MPI_Wtime()-time00 + call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + time_reduce=time_reduce+MPI_Wtime()-time00 +#define DEBUG +#ifdef DEBUG + write (iout,*) "gloc_sc after reduce" + do i=1,nres + do j=1,1 + write (iout,*) i,j,gloc_sc(j,i,icg) + enddo + enddo +#endif +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc after reduce" do i=1,4*nres @@ -5706,29 +5734,52 @@ c amino-acid residues. 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 - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1) cycle + do i=itau_start,itau_end esccor_ii=0.0D0 - itori=iabs(itype(i-2)) - itori1=iabs(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) + & 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) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci + enddo !intertyp enddo + return end c---------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index effd955..2c670f2 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -345,6 +345,7 @@ C------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.VAR' include 'COMMON.MD' + include 'COMMON.SCCOR' C C Initialize Cartesian-coordinate gradient C @@ -380,6 +381,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_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index aa5dfb9..ed4fb21 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -591,6 +591,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 @@ -1119,6 +1122,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 diff --git a/source/unres/src_MD-M/int_to_cart.f b/source/unres/src_MD-M/int_to_cart.f index 9208a80..f413622 100644 --- a/source/unres/src_MD-M/int_to_cart.f +++ b/source/unres/src_MD-M/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) @@ -113,6 +113,160 @@ c The side-chain vector derivatives enddo endif 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_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F index 33d4a50..369a4f0 100644 --- a/source/unres/src_MD-M/intcartderiv.F +++ b/source/unres/src_MD-M/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), @@ -53,7 +54,42 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i) if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint 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 @@ -123,6 +159,228 @@ 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_MD-M/kinetic_lesyng.f b/source/unres/src_MD-M/kinetic_lesyng.f index 44ea85b..db959b3 100644 --- a/source/unres/src_MD-M/kinetic_lesyng.f +++ b/source/unres/src_MD-M/kinetic_lesyng.f @@ -44,7 +44,7 @@ c to the velocities of the first Calpha. incr(j)=d_t(j,0) enddo do i=nnt,nct - iti=itype(i) + iti=iabs(itype(i)) if (itype(i).eq.10) then do j=1,3 v(j)=incr(j) diff --git a/source/unres/src_MD-M/lagrangian_lesyng.F b/source/unres/src_MD-M/lagrangian_lesyng.F index 9611f90..f8834ea 100644 --- a/source/unres/src_MD-M/lagrangian_lesyng.F +++ b/source/unres/src_MD-M/lagrangian_lesyng.F @@ -217,12 +217,12 @@ c Diagonal elements of the dX part of A and the respective friction coefficient ind=ind+1 ii = ind+m iti=itype(i) - massvec(ii)=msc(iti) + massvec(ii)=msc(iabs(iti)) if (iti.ne.10 .and. iti.ne.ntyp1) then ind1=ind1+1 ii1= ind1+m1 A(ii,ii1)=1.0d0 - Gmat(ii1,ii1)=ISC(iti) + Gmat(ii1,ii1)=ISC(iabs(iti)) endif enddo c Off-diagonal elements of the dX part of A diff --git a/source/unres/src_MD-M/moments.f b/source/unres/src_MD-M/moments.f index 5992773..983ce36 100644 --- a/source/unres/src_MD-M/moments.f +++ b/source/unres/src_MD-M/moments.f @@ -43,10 +43,10 @@ c calculating the center of the mass of the protein M_SC=0.0d0 do i=nnt,nct iti=iabs(itype(i)) - M_SC=M_SC+msc(iti) + M_SC=M_SC+msc(iabs(iti)) inres=i+nres do j=1,3 - cm(j)=cm(j)+msc(iti)*c(j,inres) + cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) enddo enddo do j=1,3 @@ -71,12 +71,12 @@ c calculating the center of the mass of the protein do j=1,3 pr(j)=c(j,inres)-cm(j) enddo - Im(1,1)=Im(1,1)+msc(iti)*(pr(2)*pr(2)+pr(3)*pr(3)) - Im(1,2)=Im(1,2)-msc(iti)*pr(1)*pr(2) - Im(1,3)=Im(1,3)-msc(iti)*pr(1)*pr(3) - Im(2,3)=Im(2,3)-msc(iti)*pr(2)*pr(3) - Im(2,2)=Im(2,2)+msc(iti)*(pr(3)*pr(3)+pr(1)*pr(1)) - Im(3,3)=Im(3,3)+msc(iti)*(pr(1)*pr(1)+pr(2)*pr(2)) + Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2) + Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3) + Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3) + Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1)) + Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct-1 @@ -262,7 +262,7 @@ c Calculate the angular momentum c write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), c & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) do j=1,3 - L(j)=L(j)+msc(iti)*vp(j) + L(j)=L(j)+msc(iabs(iti))*vp(j) enddo c write (iout,*) "L",(l(j),j=1,3) if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 9e1614b..33f017c 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -600,32 +600,129 @@ C Martix of D parameters for two dimesional fourier series enddo endif #endif +C Read of Side-chain backbone correlation parameters +C Modified 11 May 2012 by Adasko +CCC 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,ntyp - do j=1,ntyp - 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) + do i=-ntyp,-1 + isccortyp(i)=-isccortyp(-i) + 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 + nterm_sccor(-i,j)=nterm_sccor(i,j) + nterm_sccor(-i,-j)=nterm_sccor(i,j) + 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) + if (j.eq.iscprol) then + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 + & +v2sccor(k,l,i,j)*dsqrt(0.75d0) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 + & +v1sccor(k,l,i,j)*dsqrt(0.75d0) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + else + if (i.eq.isccortyp(10)) then + v1sccor(k,l,i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + else + if (j.eq.isccortyp(10)) then + v1sccor(k,l,-i,j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j) + else + v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j) + v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j) + v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j) + v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j) + v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j) + endif + endif + 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,ntyp - do j=1,ntyp + 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) + 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 endif + C C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local C interaction energy of the Gly, Ala, and Pro prototypes. diff --git a/source/unres/src_MD-M/stochfric.F b/source/unres/src_MD-M/stochfric.F index 6427551..cda93c3 100644 --- a/source/unres/src_MD-M/stochfric.F +++ b/source/unres/src_MD-M/stochfric.F @@ -39,7 +39,7 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i).ne.10) then do j=1,3 d_t_work(ind+j)=d_t(j,i+nres) enddo @@ -68,7 +68,7 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i).ne.10) then do j=1,3 friction(j,i+nres)=fric_work(ind+j) enddo @@ -185,7 +185,11 @@ c----------------------------------------------------- enddo x=0.0d0 +#ifdef MPI time00=MPI_Wtime() +#else + time00=tcpu() +#endif c Compute the stochastic forces acting on bodies. Store in force. do i=nnt,nct-1 sig=stdforcp(i) @@ -203,7 +207,11 @@ c Compute the stochastic forces acting on bodies. Store in force. force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2) enddo enddo +#ifdef MPI time_fsample=time_fsample+MPI_Wtime()-time00 +#else + time_fsample=time_fsample+tcpu()-time00 +#endif c Compute the stochastic forces acting on virtual-bond vectors. do j=1,3 ff(j)=0.0d0 @@ -226,7 +234,7 @@ c Compute the stochastic forces acting on virtual-bond vectors. stochforc(j,0)=ff(j)+force(j,nnt+nres) enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i).ne.10) then do j=1,3 stochforc(j,i+nres)=force(j,i+nres) enddo @@ -244,7 +252,7 @@ c Compute the stochastic forces acting on virtual-bond vectors. ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i).ne.10) then do j=1,3 stochforcvec(ind+j)=stochforc(j,i+nres) enddo @@ -309,7 +317,9 @@ c Compute the stochastic forces acting on virtual-bond vectors. c------------------------------------------------------------------ subroutine setup_fricmat implicit real*8 (a-h,o-z) +#ifdef MPI include 'mpif.h' +#endif include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -355,7 +365,6 @@ c Load the friction coefficients corresponding to peptide groups c Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 - gamsc(ntyp1)=1.0d0 do i=nnt,nct ind=ind+1 ii = ind+m @@ -478,22 +487,39 @@ c enddo if (nfgtasks.gt.1) then if (fg_rank.eq.0) then c The matching BROADCAST for fg processors is called in ERGASTULUM +#ifdef MPI time00=MPI_Wtime() +#else + time00=tcpu() +#endif call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR) +#ifdef MPI time_Bcast=time_Bcast+MPI_Wtime()-time00 +#else + time_Bcast=time_Bcast+tcpu()-time00 +#endif c print *,"Processor",myrank, c & " BROADCAST iorder in SETUP_FRICMAT" endif c licznik=licznik+1 c write (iout,*) "setup_fricmat licznik",licznik +#ifdef MPI time00=MPI_Wtime() +#else + time00=tcpu() +#endif c Scatter the friction matrix call MPI_Scatterv(fricmat(1,1),nginv_counts(0), & nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1), & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) - time_scatter=time_scatter+MPI_Wtime()-time00 #ifdef TIMING +#ifdef MPI + time_scatter=time_scatter+MPI_Wtime()-time00 time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00 +#else + time_scatter=time_scatter+tcpu()-time00 + time_scatter_fmat=time_scatter_fmat+tcpu()-time00 +#endif #endif do i=1,dimen do j=1,2*my_ng_count @@ -532,7 +558,7 @@ c include 'COMMON.NAMES' double precision radius(maxres2),gamvec(maxres6) parameter (twosix=1.122462048309372981d0) - logical lprn /.true./ + logical lprn /.false./ c c determine new friction coefficients every few SD steps c diff --git a/source/unres/src_MD/COMMON.SCCOR b/source/unres/src_MD/COMMON.SCCOR index 6e39f93..e29cb4c 100644 --- a/source/unres/src_MD/COMMON.SCCOR +++ b/source/unres/src_MD/COMMON.SCCOR @@ -4,9 +4,9 @@ cc Parameters of the SCCOR term & dcostau,dsintau,dtauangle,dcosomicron, & domicron integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor - common/sccor/v1sccor(maxterm_sccor,3,20,20), - & v2sccor(maxterm_sccor,3,20,20), - & v0sccor(maxterm_sccor,20), + 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),