From: Adam Liwo Date: Sun, 20 Dec 2015 19:02:54 +0000 (+0100) Subject: Merge branch 'prerelease-3.2.1' of mmka.chem.univ.gda.pl:unres into prerelease-3.2.1 X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=924b9537a57bfdc673ec7bbb0f978024201b6021;hp=95abe93df3a12f9cf8d94b2ed69c14cc4b3363e0;p=unres.git Merge branch 'prerelease-3.2.1' of mmka.chem.univ.gda.pl:unres into prerelease-3.2.1 --- diff --git a/source/unres/src_CSA/CMakeLists.txt b/source/unres/src_CSA/CMakeLists.txt index 340fd7c..2f035cb 100644 --- a/source/unres/src_CSA/CMakeLists.txt +++ b/source/unres/src_CSA/CMakeLists.txt @@ -64,6 +64,7 @@ set(UNRES_CSA_SRC0 TMscore_subroutine.f together.F unres_csa.F + gnmr1.f ) set(UNRES_CSA_SRC3 energy_p_new_barrier.F gradient_p.F ) diff --git a/source/unres/src_CSA/COMMON.SBRIDGE b/source/unres/src_CSA/COMMON.SBRIDGE index 4cc80c8..953b380 100644 --- a/source/unres/src_CSA/COMMON.SBRIDGE +++ b/source/unres/src_CSA/COMMON.SBRIDGE @@ -4,8 +4,8 @@ & ns,nss,nfree,iss(maxss) double precision dhpb,forcon integer ihpb,jhpb,nhpb - common /links/ dhpb(maxdim),forcon(maxdim),ihpb(maxdim), - & jhpb(maxdim),nhpb + common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), + & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb double precision weidis common /restraints/ weidis integer link_start,link_end diff --git a/source/unres/src_CSA/Makefile_MPICH_bluegene b/source/unres/src_CSA/Makefile_MPICH_bluegene index bcfd204..72b069a 100644 --- a/source/unres/src_CSA/Makefile_MPICH_bluegene +++ b/source/unres/src_CSA/Makefile_MPICH_bluegene @@ -48,7 +48,7 @@ object = unres_csa.o arcos.o cartprint.o chainbuild.o initialize_p.o \ together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \ indexx.o prng_32.o contact.o gen_rand_conf.o \ sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \ - distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o + distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o no_option: diff --git a/source/unres/src_CSA/Makefile_MPICH_gfortran b/source/unres/src_CSA/Makefile_MPICH_gfortran index 529d775..8132171 100644 --- a/source/unres/src_CSA/Makefile_MPICH_gfortran +++ b/source/unres/src_CSA/Makefile_MPICH_gfortran @@ -37,7 +37,7 @@ object = unres_csa.o arcos.o cartprint.o chainbuild.o initialize_p.o \ together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \ indexx.o prng_32.o contact.o gen_rand_conf.o \ sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \ - distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o + distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o no_option: diff --git a/source/unres/src_CSA/Makefile_MPICH_ifort b/source/unres/src_CSA/Makefile_MPICH_ifort index 8a8aa2b..dfa0875 100644 --- a/source/unres/src_CSA/Makefile_MPICH_ifort +++ b/source/unres/src_CSA/Makefile_MPICH_ifort @@ -38,7 +38,7 @@ object = unres_csa.o arcos.o cartprint.o chainbuild.o initialize_p.o \ together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \ indexx.o prng_32.o contact.o gen_rand_conf.o \ sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \ - distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o + distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o no_option: diff --git a/source/unres/src_CSA/energy_p_new_barrier.F b/source/unres/src_CSA/energy_p_new_barrier.F index 4ca78f8..2121254 100644 --- a/source/unres/src_CSA/energy_p_new_barrier.F +++ b/source/unres/src_CSA/energy_p_new_barrier.F @@ -4239,26 +4239,72 @@ C iii and jjj point to the residues for which the distance is assigned. iii=ii jjj=jj endif -cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj +c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, +c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds + if (i.le.nss) then + if (ii.gt.nres + & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + dd=dist(ii,jj) + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "beta reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + do k=1,3 + ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) + ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) + enddo else C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) C Calculate the contribution to energy. ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis C C Evaluate gradient. C fac=waga*rdis/dd + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 diff --git a/source/unres/src_CSA/gnmr1.f b/source/unres/src_CSA/gnmr1.f new file mode 100644 index 0000000..905e746 --- /dev/null +++ b/source/unres/src_CSA/gnmr1.f @@ -0,0 +1,43 @@ + double precision function gnmr1(y,ymin,ymax) + implicit none + double precision y,ymin,ymax + double precision wykl /4.0d0/ + if (y.lt.ymin) then + gnmr1=(ymin-y)**wykl/wykl + else if (y.gt.ymax) then + gnmr1=(y-ymax)**wykl/wykl + else + gnmr1=0.0d0 + endif + return + end +c------------------------------------------------------------------------------ + double precision function gnmr1prim(y,ymin,ymax) + implicit none + double precision y,ymin,ymax + double precision wykl /4.0d0/ + if (y.lt.ymin) then + gnmr1prim=-(ymin-y)**(wykl-1) + else if (y.gt.ymax) then + gnmr1prim=(y-ymax)**(wykl-1) + else + gnmr1prim=0.0d0 + endif + return + end +c------------------------------------------------------------------------------ + double precision function harmonic(y,ymax) + implicit none + double precision y,ymax + double precision wykl /2.0d0/ + harmonic=(y-ymax)**wykl + return + end +c------------------------------------------------------------------------------- + double precision function harmonicprim(y,ymax) + double precision y,ymin,ymax + double precision wykl /2.0d0/ + harmonicprim=(y-ymax)*wykl + return + end +c--------------------------------------------------------------------------------- diff --git a/source/unres/src_CSA/readrtns_csa.F b/source/unres/src_CSA/readrtns_csa.F index e88a67e..af05ce6 100644 --- a/source/unres/src_CSA/readrtns_csa.F +++ b/source/unres/src_CSA/readrtns_csa.F @@ -665,11 +665,6 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup - call flush(iout) - if (constr_dist.gt.0) call read_dist_constr -c write (iout,*) "After read_dist_constr nhpb",nhpb - call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co if (pdbref) then @@ -686,6 +681,13 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb enddo endif endif +c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + if (constr_dist.gt.0) then + call read_dist_constr + endif + if (nhpb.gt.0) call hpb_partition +c write (iout,*) "After read_dist_constr nhpb",nhpb +c call flush(iout) if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@ -1668,9 +1670,9 @@ c---------------------------------------------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard - write (iout,*) "Calling read_dist_constr" - write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup - call flush(iout) +c write (iout,*) "Calling read_dist_constr" +c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +c call flush(iout) call card_concat(controlcard) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) @@ -1689,6 +1691,18 @@ c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo + if (.not.refstr .and. nfrag.gt.0) then + write (iout,*) + & "ERROR: no reference structure to compute distance restraints" + write (iout,*) + & "Restraints must be specified explicitly (NDIST=number)" + stop + endif + if (nfrag.lt.2 .and. npair.gt.0) then + write (iout,*) "ERROR: Less than 2 fragments specified", + & " but distance restraints between pairs requested" + stop + endif call flush(iout) do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup @@ -1699,7 +1713,7 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) if (wfrag_(i).gt.0.0d0) then do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) - write (iout,*) "j",j," k",k +c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=nhpb+1 @@ -1763,21 +1777,29 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) endif enddo do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif + enddo #ifdef MPI - if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) -#else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + if (.not.out1file .or. me.eq.king) then #endif - endif + do i=1,nhpb + write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", + & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) enddo call flush(iout) +#ifdef MPI + endif +#endif return end c------------------------------------------------------------------------------- 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 539102b..eb1cbd4 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -718,7 +718,7 @@ c enddo do i=1,4*nres glocbuf(i)=gloc(i,icg) enddo -#define DEBUG +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc_sc before reduce" do i=1,nres @@ -747,7 +747,7 @@ c enddo 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 +#undef DEBUG #ifdef DEBUG write (iout,*) "gloc_sc after reduce" do i=1,nres @@ -5648,6 +5648,7 @@ C 6/23/01 Compute double torsional energy include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging lprn=.false. @@ -5658,6 +5659,7 @@ C write(iout,*) "a tu??" if (itype(i-2).eq.21 .or. itype(i-1).eq.21 & .or. itype(i).eq.21 .or. itype(i+1).eq.21 & .or. itype(i-3).eq.ntyp1) cycle + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5677,6 +5679,8 @@ C Regular cosine and sine terms sinphi2=dsin(j*phii1) etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ & v2cij*cosphi2+v2sij*sinphi2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo @@ -5692,12 +5696,17 @@ C Regular cosine and sine terms sinphi1m2=dsin(l*phii-(k-l)*phii1) etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'etor_d',i,etors_d_ii gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 enddo diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F index 47583a4..9c25867 100644 --- a/source/unres/src_MD/energy_p_new_barrier.F +++ b/source/unres/src_MD/energy_p_new_barrier.F @@ -4493,6 +4493,8 @@ c do i=ibondp_start,ibondp_end diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + if (energy_dec) write (iout,'(a7,i5,4f7.3)') + & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) @@ -4511,6 +4513,9 @@ c diff=vbld(i+nres)-vbldsc0(1,iti) c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, c & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (energy_dec) write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -4946,6 +4951,8 @@ C if (i.gt.3) then & 'ebe', i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 gloc(nphi+i-2,icg)=wang*dethetai @@ -5887,12 +5894,14 @@ C 6/23/01 Compute double torsional energy include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 do i=iphid_start,iphid_end + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5911,6 +5920,8 @@ c lprn=.true. sinphi2=dsin(j*phii1) etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ & v2cij*cosphi2+v2sij*sinphi2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo @@ -5926,12 +5937,17 @@ c lprn=.true. sinphi1m2=dsin(l*phii-(k-l)*phii1) etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'etor_d',i,etors_d_ii gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 c write (iout,*) "gloci", gloc(i-3,icg)