From 56b6e687a929b4935fe39b4abd8beb4e05265199 Mon Sep 17 00:00:00 2001 From: Adam Kazimierz Sieradzan Date: Thu, 27 Jun 2013 21:00:10 -0400 Subject: [PATCH] Bugfix for SS wham and introduction SS to cluster analysis --- source/cluster/wham/src/CMakeLists.txt | 4 +- source/cluster/wham/src/energy_p_new.F | 36 +- .../cluster/wham/src/include_unres/COMMON.SBRIDGE | 19 +- source/cluster/wham/src/initialize_p.F | 4 + source/cluster/wham/src/main_clust.F | 1 - source/cluster/wham/src/parmread.F | 23 +- source/cluster/wham/src/probabl.F | 31 +- source/cluster/wham/src/read_coords.F | 53 +- source/cluster/wham/src/readrtns.F | 51 + source/cluster/wham/src/ssMD.F | 1961 ++++++++++++++++++++ source/wham/src/COMMON.ALLPARM | 2 + source/wham/src/bxread.F | 2 + source/wham/src/cxread.F | 13 +- source/wham/src/enecalc1.F | 66 +- source/wham/src/energy_p_new.F | 2 +- source/wham/src/geomout.F | 4 + source/wham/src/make_ensemble1.F | 14 +- source/wham/src/parmread.F | 32 +- source/wham/src/store_parm.F | 4 + source/wham/src/wham_calc1.F | 12 +- 20 files changed, 2253 insertions(+), 81 deletions(-) create mode 100644 source/cluster/wham/src/ssMD.F diff --git a/source/cluster/wham/src/CMakeLists.txt b/source/cluster/wham/src/CMakeLists.txt index 0851e71..760269e 100644 --- a/source/cluster/wham/src/CMakeLists.txt +++ b/source/cluster/wham/src/CMakeLists.txt @@ -35,6 +35,7 @@ set(UNRES_CLUSTER_WHAM_SRC0 rescode.f setup_var.f srtclust.f + ssMD.F timing.F track.F wrtclust.f @@ -50,6 +51,7 @@ set(UNRES_CLUSTER_WHAM_PP_SRC probabl.F read_coords.F readrtns.F + ssMD.F timing.F track.F work_partition.F @@ -90,7 +92,7 @@ endif(UNRES_MD_FF STREQUAL "GAB") #========================================= # Additional flags #========================================= -set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN" ) +set(CPPFLAGS "${CPPFLAGS} -DUNRES -DISNAN -DCLUST" ) #========================================= # Compiler specific flags diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index fd30820..636f983 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -107,7 +107,7 @@ C #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*fact(1)*ees+wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d @@ -115,7 +115,7 @@ C #else etot=wsc*evdw+wscp*evdw2+welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d @@ -152,6 +152,7 @@ C energia(18)=estr energia(19)=esccor energia(20)=edihcnstr +cc if (dyn_ss) call dyn_set_nss c detecting NaNQ i=0 #ifdef WINPGI @@ -723,6 +724,7 @@ c include "DIMENSIONS.COMPAR" include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant @@ -748,6 +750,12 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij +c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') +c & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=itype(j) dscj_inv=vbld_inv(j+nres) @@ -830,6 +838,7 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad endif + ENDIF ! SSBOND enddo ! j enddo ! iint enddo ! i @@ -854,6 +863,7 @@ c include "DIMENSIONS.COMPAR" include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' common /srutu/ icall logical lprn integer icant @@ -879,6 +889,13 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) +C in case of diagnostics write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij +c if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') +c & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=itype(j) dscj_inv=vbld_inv(j+nres) @@ -961,6 +978,7 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -2822,10 +2840,13 @@ 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 (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check 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 cd write (iout,*) "eij",eij + endif else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction dd=dist(ii,jj) @@ -2957,7 +2978,7 @@ C deltat12=om2-om1+2.0d0 cosphi=om12-om1*om2 eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) - & +akct*deltad*deltat12 + & +akct*deltad*deltat12+ebr & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi c write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, c & " akct",akct," deltad",deltad," deltat",deltat1,deltat2, @@ -4869,6 +4890,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding C Set lprn=.true. for debugging lprn=.false. eturn6=0.0d0 + ecorr6=0.0d0 #ifdef MPL n_corr=0 n_corr1=0 @@ -5045,10 +5067,10 @@ cd write(2,*)'wcorr6',wcorr6,' wturn6',wturn6 cd write(2,*)'ijkl',i,j,i+1,j1 if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3 & .or. wturn6.eq.0.0d0))then -cd write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 - ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk) -cd write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, -cd & 'ecorr6=',ecorr6 +c write (iout,*) '******ecorr6: i,j,i+1,j1',i,j,i+1,j1 +c ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk) +c write (iout,*) 'ecorr',ecorr,' ecorr5=',ecorr5, +c & 'ecorr6=',ecorr6, wcorr6 cd write (iout,'(4e15.5)') sred_geom, cd & dabs(eello4(i,j,i+1,j1,jj,kk)), cd & dabs(eello5(i,j,i+1,j1,jj,kk)), diff --git a/source/cluster/wham/src/include_unres/COMMON.SBRIDGE b/source/cluster/wham/src/include_unres/COMMON.SBRIDGE index 7bba010..f866aa7 100644 --- a/source/cluster/wham/src/include_unres/COMMON.SBRIDGE +++ b/source/cluster/wham/src/include_unres/COMMON.SBRIDGE @@ -1,10 +1,17 @@ - double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,dhpb, - & dhpb1,forcon,weidis - integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end, - & ibecarb - common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,ns,nss, - & nfree,iss(maxss) + double precision ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss + integer ns,nss,nfree,iss + common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, + & ns,nss,nfree,iss(maxss) + double precision dhpb,dhpb1,forcon + integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb 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 common /links_split/ link_start,link_end + double precision Ht,dyn_ssbond_ij + logical dyn_ss,dyn_ss_mask + common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres), + & idssb(maxdim),jdssb(maxdim), + & Ht,dyn_ss,dyn_ss_mask(maxres) diff --git a/source/cluster/wham/src/initialize_p.F b/source/cluster/wham/src/initialize_p.F index 37e0bf9..f8b9426 100644 --- a/source/cluster/wham/src/initialize_p.F +++ b/source/cluster/wham/src/initialize_p.F @@ -155,6 +155,9 @@ C Initialize the bridge arrays ihpb(i)=0 jhpb(i)=0 enddo + do i=1,maxres + dyn_ss_mask(i)=.false. + enddo C C Initialize timing. C @@ -291,6 +294,7 @@ cd & (ihpb(i),jhpb(i),i=1,nss) do ii=1,nss if (ihpb(ii).eq.i+nres) then scheck=.true. + if (dyn_ss) go to 10 jj=jhpb(ii)-nres goto 10 endif diff --git a/source/cluster/wham/src/main_clust.F b/source/cluster/wham/src/main_clust.F index 4f50091..a0ae38f 100644 --- a/source/cluster/wham/src/main_clust.F +++ b/source/cluster/wham/src/main_clust.F @@ -107,7 +107,6 @@ c if (refstr) call read_ref_structure(*30) #ifdef MPI call work_partition(.true.,ncon) #endif - call probabl(iT,ncon_work,ncon,*20) if (ncon_work.lt.2) then diff --git a/source/cluster/wham/src/parmread.F b/source/cluster/wham/src/parmread.F index 7f6a145..656d219 100644 --- a/source/cluster/wham/src/parmread.F +++ b/source/cluster/wham/src/parmread.F @@ -810,7 +810,7 @@ C C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +C ebr=-5.50D0 c c Old arbitrary potential - commented out. c @@ -821,19 +821,12 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G** c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 - - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, - & ' v3ss:',v3ss +C D0CM = 3.78d0 +C AKCM = 15.1d0 +C AKTH = 11.0d0 +C AKCT = 12.0d0 +C V1SS =-1.08d0 +C V2SS = 7.61d0 +C V3SS = 13.7d0 return end diff --git a/source/cluster/wham/src/probabl.F b/source/cluster/wham/src/probabl.F index 4c49d99..7fcd29b 100644 --- a/source/cluster/wham/src/probabl.F +++ b/source/cluster/wham/src/probabl.F @@ -36,11 +36,20 @@ c do i=1,ncon c write (iout,*) i,list_conf(i) c enddo +c do i=1,ncon +c write(iout,*) "entrop before", entfac(i),i +c enddo + #ifdef MPI write (iout,*) me," indstart",indstart(me)," indend",indend(me) call daread_ccoords(indstart(me),indend(me)) #endif +c do i=1,ncon +c write(iout,*) "entrop after", entfac(i),i +c enddo + c write (iout,*) "ncon",ncon + temper=1.0d0/(beta_h(ib)*1.987D-3) c write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper c quot=1.0d0/(T0*beta_h(ib)*1.987D-3) @@ -102,11 +111,13 @@ c write (iout,*) "i",i," ii",ii call int_from_cart1(.false.) call etotal(energia(0),fT) totfree(i)=energia(0) +c#define DEBUG #ifdef DEBUG - write (iout,*) i," energia",(energia(j),j=0,21) + write (iout,*) i," energia",(energia(j),j=0,20) call enerprint(energia(0),ft) call flush(iout) #endif +c#undef DEBUG do k=1,max_ene enetb(k,i)=energia(k) enddo @@ -129,6 +140,7 @@ c write (iout,*) "i",i," ii",ii ecorr=enetb(4,i) ecorr5=enetb(5,i) ecorr6=enetb(6,i) +cc if (wcorr6.eq.0) ecorr6=0.0d0 eel_loc=enetb(7,i) eello_turn3=enetb(8,i) eello_turn4=enetb(9,i) @@ -144,7 +156,7 @@ c write (iout,*) "i",i," ii",ii #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+ft(1)*welec*ees+wvdwpp*evdw1 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc @@ -153,7 +165,7 @@ c write (iout,*) "i",i," ii",ii #else etot=wsc*evdw+wscp*evdw2+ft(1)*welec*(ees+evdw1) & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr @@ -163,7 +175,8 @@ c write (iout,*) "i",i," ii",ii Fdimless(i)=beta_h(ib)*etot+entfac(ii) totfree(i)=etot #ifdef DEBUG - write (iout,*) i,ii,ib, + + write (iout,*) "etrop", i,ii,ib, & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), & entfac(ii),Fdimless(i) #endif @@ -183,6 +196,7 @@ c write (iout,*) "i",i," ii",ii & MPI_COMM_WORLD, IERROR) if (me.eq.Master) then #endif +c#define DEBUG #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ncon @@ -196,24 +210,27 @@ c write (iout,*) "i",i," ii",ii write (iout,*) i,list_conf(i),fdimless(i) enddo #endif +c#undef DEBUG do i=1,ncon totfree(i)=fdimless(i) enddo qfree=0.0d0 do i=1,ncon - qfree=qfree+exp(-fdimless(i)+fdimless(1)) + qfree=qfree+dexp(dble(-fdimless(i)+fdimless(1))) enddo c write (iout,*) "qfree",qfree nlist=1 sumprob=0.0 do i=1,min0(ncon,maxstr_proc)-1 - sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree + sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree +c#define DEBUG #ifdef DEBUG - write (iout,*) i,ib,beta_h(ib), + write (iout,*) 'i=',i,ib,beta_h(ib), & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i), & totfree(list_conf(i)), & -entfac(list_conf(i)),fdimless(i),sumprob #endif +c#undef DEBUG if (sumprob.gt.prob_limit) goto 122 c if (sumprob.gt.1.00d0) goto 122 nlist=nlist+1 diff --git a/source/cluster/wham/src/read_coords.F b/source/cluster/wham/src/read_coords.F index 2a21cbe..cf98db7 100644 --- a/source/cluster/wham/src/read_coords.F +++ b/source/cluster/wham/src/read_coords.F @@ -218,10 +218,19 @@ c call flush(iout) call xdrfint_(ixdrf, nss, iret) if (iret.eq.0) goto 101 do j=1,nss - call xdrfint_(ixdrf, ihpb(j), iret) - if (iret.eq.0) goto 101 - call xdrfint_(ixdrf, jhpb(j), iret) - if (iret.eq.0) goto 101 +cc if (dyn_ss) then +cc call xdrfint_(ixdrf, idssb(j), iret) +cc if (iret.eq.0) goto 101 +cc call xdrfint_(ixdrf, jdssb(j), iret) +cc if (iret.eq.0) goto 101 +cc idssb(j)=idssb(j)-nres +cc jdssb(j)=jdssb(j)-nres +cc else + call xdrfint_(ixdrf, ihpb(j), iret) + if (iret.eq.0) goto 101 + call xdrfint_(ixdrf, jhpb(j), iret) + if (iret.eq.0) goto 101 +cc endif enddo call xdrffloat_(ixdrf,reini,iret) if (iret.eq.0) goto 101 @@ -243,10 +252,20 @@ c write (iout,*) "nss",nss call flush(iout) if (iret.eq.0) goto 101 do k=1,nss - call xdrfint(ixdrf, ihpb(k), iret) - if (iret.eq.0) goto 101 - call xdrfint(ixdrf, jhpb(k), iret) - if (iret.eq.0) goto 101 +cc if (dyn_ss) then +cc call xdrfint(ixdrf, idssb(k), iret) +cc if (iret.eq.0) goto 101 +cc call xdrfint(ixdrf, jdssb(k), iret) +cc if (iret.eq.0) goto 101 +cc idssb(k)=idssb(k)-nres +cc jdssb(k)=jdssb(k)-nres +cc write(iout,*) "TUTU", idssb(k),jdssb(k) +cc else + call xdrfint(ixdrf, ihpb(k), iret) + if (iret.eq.0) goto 101 + call xdrfint(ixdrf, jhpb(k), iret) + if (iret.eq.0) goto 101 +cc endif enddo call xdrffloat(ixdrf,reini,iret) if (iret.eq.0) goto 101 @@ -258,7 +277,9 @@ c write (iout,*) "nss",nss if (iret.eq.0) goto 101 #endif energy(jj+1)=reini - entfac(jj+1)=refree +cc write(iout,*) 'reini=', reini, jj+1 + entfac(jj+1)=dble(refree) +cc write(iout,*) 'refree=', refree,jj+1 rmstb(jj+1)=rmsdev do k=1,nres do l=1,3 @@ -639,10 +660,17 @@ c write (iout,*) "Reading binary file, record",iii," ii",ii call flush(iout) #endif + if (dyn_ss) then + read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), + & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), +c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss), + & entfac(ii),rmstb(ii) + else read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss), & entfac(ii),rmstb(ii) + endif #ifdef DEBUG write (iout,*) ii,iii,ij,entfac(ii) write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) @@ -696,10 +724,17 @@ c write (iout,*) "Writing binary file, record",iii," ii",ii call flush(iout) #endif + if (dyn_ss) then + write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), + & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), +c & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)) + & entfac(ii),rmstb(ii) + else write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres), & nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)), & entfac(ii),rmstb(ii) + endif #ifdef DEBUG write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres, diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index 8e63ff8..c40fcbb 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -30,6 +30,7 @@ C refstr=(index(controlcard,'REFSTR').gt.0) write (iout,*) "REFSTR",refstr pdbref=(index(controlcard,'PDBREF').gt.0) + dyn_ss=(index(controlcard,'DYN_SS').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) min_var=(index(controlcard,'MINVAR').gt.0) @@ -126,9 +127,43 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif + write (iout,'(/a)') "Disulfide bridge parameters:" + write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr + write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth + write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm + write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct + write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, + & ' v3ss:',v3ss + write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp @@ -431,6 +466,22 @@ c forcon(i)=fbr enddo endif endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. +c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU" + enddo + endif + print *, "Leaving brigde read" return end c---------------------------------------------------------------------------- diff --git a/source/cluster/wham/src/ssMD.F b/source/cluster/wham/src/ssMD.F new file mode 100644 index 0000000..70ed6fd --- /dev/null +++ b/source/cluster/wham/src/ssMD.F @@ -0,0 +1,1961 @@ +c---------------------------------------------------------------------------- + subroutine check_energies +c implicit none + +c Includes + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + include 'COMMON.LOCAL' + include 'COMMON.GEO' + +c External functions + double precision ran_number + external ran_number + +c Local variables + integer i,j,k,l,lmax,p,pmax + double precision rmin,rmax + double precision eij + + double precision d + double precision wi,rij,tj,pj + + +c return + + i=5 + j=14 + + d=dsc(1) + rmin=2.0D0 + rmax=12.0D0 + + lmax=10000 + pmax=1 + + do k=1,3 + c(k,i)=0.0D0 + c(k,j)=0.0D0 + c(k,nres+i)=0.0D0 + c(k,nres+j)=0.0D0 + enddo + + do l=1,lmax + +ct wi=ran_number(0.0D0,pi) +c wi=ran_number(0.0D0,pi/6.0D0) +c wi=0.0D0 +ct tj=ran_number(0.0D0,pi) +ct pj=ran_number(0.0D0,pi) +c pj=ran_number(0.0D0,pi/6.0D0) +c pj=0.0D0 + + do p=1,pmax +ct rij=ran_number(rmin,rmax) + + c(1,j)=d*sin(pj)*cos(tj) + c(2,j)=d*sin(pj)*sin(tj) + c(3,j)=d*cos(pj) + + c(3,nres+i)=-rij + + c(1,i)=d*sin(wi) + c(3,i)=-rij-d*cos(wi) + + do k=1,3 + dc(k,nres+i)=c(k,nres+i)-c(k,i) + dc_norm(k,nres+i)=dc(k,nres+i)/d + dc(k,nres+j)=c(k,nres+j)-c(k,j) + dc_norm(k,nres+j)=dc(k,nres+j)/d + enddo + + call dyn_ssbond_ene(i,j,eij) + enddo + enddo + + call exit(1) + + return + end + +C----------------------------------------------------------------------------- + + subroutine dyn_ssbond_ene(resi,resj,eij) +c implicit none + +c Includes + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' +#endif +#endif + +c External functions + double precision h_base + external h_base + +c Input arguments + integer resi,resj + +c Output arguments + double precision eij + +c Local variables + logical havebond +c integer itypi,itypj,k,l + double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi + double precision sig0ij,ljd,sig,fac,e1,e2 + double precision dcosom1(3),dcosom2(3),ed + double precision pom1,pom2 + double precision ljA,ljB,ljXs + double precision d_ljB(1:3) + double precision ssA,ssB,ssC,ssXs + double precision ssxm,ljxm,ssm,ljm + double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3) + double precision f1,f2,h1,h2,hd1,hd2 + double precision omega,delta_inv,deltasq_inv,fac1,fac2 +c-------FIRST METHOD + double precision xm,d_xm(1:3) +c-------END FIRST METHOD +c-------SECOND METHOD +c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3) +c-------END SECOND METHOD + +c-------TESTING CODE + logical checkstop,transgrad + common /sschecks/ checkstop,transgrad + + integer icheck,nicheck,jcheck,njcheck + double precision echeck(-1:1),deps,ssx0,ljx0 +c-------END TESTING CODE + + + i=resi + j=resj + + itypi=itype(i) + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=vbld_inv(i+nres) + + itypj=itype(j) + xj=c(1,nres+j)-c(1,nres+i) + yj=c(2,nres+j)-c(2,nres+i) + zj=c(3,nres+j)-c(3,nres+i) + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + dscj_inv=vbld_inv(j+nres) + + chi1=chi(itypi,itypj) + chi2=chi(itypj,itypi) + chi12=chi1*chi2 + chip1=chip(itypi) + chip2=chip(itypj) + chip12=chip1*chip2 + alf1=alp(itypi) + alf2=alp(itypj) + alf12=0.5D0*(alf1+alf2) + + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse +c The following are set in sc_angular +c erij(1)=xj*rij +c erij(2)=yj*rij +c erij(3)=zj*rij +c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) +c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) +c om12=dxi*dxj+dyi*dyj+dzi*dzj + call sc_angular + rij=1.0D0/rij ! Reset this so it makes sense + + sig0ij=sigma(itypi,itypj) + sig=sig0ij*dsqrt(1.0D0/sigsq) + + ljXs=sig-sig0ij + ljA=eps1*eps2rt**2*eps3rt**2 + ljB=ljA*bb(itypi,itypj) + ljA=ljA*aa(itypi,itypj) + ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + + ssXs=d0cm + deltat1=1.0d0-om1 + deltat2=1.0d0+om2 + deltat12=om2-om1+2.0d0 + cosphi=om12-om1*om2 + ssA=akcm + ssB=akct*deltat12 + ssC=ss_depth + & +akth*(deltat1*deltat1+deltat2*deltat2) + & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi + ssxm=ssXs-0.5D0*ssB/ssA + +c-------TESTING CODE +c$$$c Some extra output +c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA +c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) +c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC +c$$$ if (ssx0.gt.0.0d0) then +c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA +c$$$ else +c$$$ ssx0=ssxm +c$$$ endif +c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ", +c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12 +c$$$ return +c-------END TESTING CODE + +c-------TESTING CODE +c Stop and plot energy and derivative as a function of distance + if (checkstop) then + ssm=ssC-0.25D0*ssB*ssB/ssA + ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + if (ssm.lt.ljm .and. + & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then + nicheck=1000 + njcheck=1 + deps=0.5d-7 + else + checkstop=.false. + endif + endif + if (.not.checkstop) then + nicheck=0 + njcheck=-1 + endif + + do icheck=0,nicheck + do jcheck=-1,njcheck + if (checkstop) rij=(ssxm-1.0d0)+ + & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps +c-------END TESTING CODE + + if (rij.gt.ljxm) then + havebond=.false. + ljd=rij-ljXs + fac=(1.0D0/ljd)**expon + e1=fac*fac*aa(itypi,itypj) + e2=fac*bb(itypi,itypj) + eij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=eij*eps3rt + eps3der=eij*eps2rt + eij=eij*eps2rt*eps3rt + + sigder=-sig/sigsq + e1=e1*eps1*eps2rt**2*eps3rt**2 + ed=-expon*(e1+eij)/ljd + sigder=ed*sigder + eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 + eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 + eom12=eij*eps1_om12+eps2der*eps2rt_om12 + & -2.0D0*alf12*eps3der+sigder*sigsq_om12 + else if (rij.lt.ssxm) then + havebond=.true. + ssd=rij-ssXs + eij=ssA*ssd*ssd+ssB*ssd+ssC + + ed=2*akcm*ssd+akct*deltat12 + pom1=akct*ssd + pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi + eom1=-2*akth*deltat1-pom1-om2*pom2 + eom2= 2*akth*deltat2+pom1-om1*pom2 + eom12=pom2 + else + omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi + + d_ssxm(1)=0.5D0*akct/ssA + d_ssxm(2)=-d_ssxm(1) + d_ssxm(3)=0.0D0 + + d_ljxm(1)=sig0ij/sqrt(sigsq**3) + d_ljxm(2)=d_ljxm(1)*sigsq_om2 + d_ljxm(3)=d_ljxm(1)*sigsq_om12 + d_ljxm(1)=d_ljxm(1)*sigsq_om1 + +c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE + xm=0.5d0*(ssxm+ljxm) + do k=1,3 + d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k)) + enddo + if (rij.lt.xm) then + havebond=.true. + ssm=ssC-0.25D0*ssB*ssB/ssA + d_ssm(1)=0.5D0*akct*ssB/ssA + d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1) + d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1) + d_ssm(3)=omega + f1=(rij-xm)/(ssxm-xm) + f2=(rij-ssxm)/(xm-ssxm) + h1=h_base(f1,hd1) + h2=h_base(f2,hd2) + eij=ssm*h1+Ht*h2 + delta_inv=1.0d0/(xm-ssxm) + deltasq_inv=delta_inv*delta_inv + fac=ssm*hd1-Ht*hd2 + fac1=deltasq_inv*fac*(xm-rij) + fac2=deltasq_inv*fac*(rij-ssxm) + ed=delta_inv*(Ht*hd2-ssm*hd1) + eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1) + eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2) + eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3) + else + havebond=.false. + ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) + d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- + + alf12/eps3rt) + d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt) + f1=(rij-ljxm)/(xm-ljxm) + f2=(rij-xm)/(ljxm-xm) + h1=h_base(f1,hd1) + h2=h_base(f2,hd2) + eij=Ht*h1+ljm*h2 + delta_inv=1.0d0/(ljxm-xm) + deltasq_inv=delta_inv*delta_inv + fac=Ht*hd1-ljm*hd2 + fac1=deltasq_inv*fac*(ljxm-rij) + fac2=deltasq_inv*fac*(rij-xm) + ed=delta_inv*(ljm*hd2-Ht*hd1) + eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1) + eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2) + eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3) + endif +c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE + +c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE +c$$$ ssd=rij-ssXs +c$$$ ljd=rij-ljXs +c$$$ fac1=rij-ljxm +c$$$ fac2=rij-ssxm +c$$$ +c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt) +c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt) +c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt) +c$$$ +c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA +c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA +c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1) +c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1) +c$$$ d_ssm(3)=omega +c$$$ +c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj) +c$$$ do k=1,3 +c$$$ d_ljm(k)=ljm*d_ljB(k) +c$$$ enddo +c$$$ ljm=ljm*ljB +c$$$ +c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC +c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB +c$$$ d_ss(2)=akct*ssd +c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega +c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega +c$$$ d_ss(3)=omega +c$$$ +c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj) +c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0) +c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1 +c$$$ do k=1,3 +c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1- +c$$$ & 2.0d0*ljB*fac1*d_ljxm(k)) +c$$$ enddo +c$$$ ljf=ljm+ljf*ljB*fac1*fac1 +c$$$ +c$$$ f1=(rij-ljxm)/(ssxm-ljxm) +c$$$ f2=(rij-ssxm)/(ljxm-ssxm) +c$$$ h1=h_base(f1,hd1) +c$$$ h2=h_base(f2,hd2) +c$$$ eij=ss*h1+ljf*h2 +c$$$ delta_inv=1.0d0/(ljxm-ssxm) +c$$$ deltasq_inv=delta_inv*delta_inv +c$$$ fac=ljf*hd2-ss*hd1 +c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac +c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1))) +c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2))) +c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3))) +c$$$ +c$$$ havebond=.false. +c$$$ if (ed.gt.0.0d0) havebond=.true. +c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE + + endif + + if (havebond) then +#ifndef CLUST +#ifndef WHAM +c if (dyn_ssbond_ij(i,j).eq.1.0d300) then +c write(iout,'(a15,f12.2,f8.1,2i5)') +c & "SSBOND_E_FORM",totT,t_bath,i,j +c endif +#endif +#endif + write(iout,*), 'DYN_SS_BOND',i,j,eij + dyn_ssbond_ij(i,j)=eij + else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then + dyn_ssbond_ij(i,j)=1.0d300 +#ifndef CLUST +#ifndef WHAM +c write(iout,'(a15,f12.2,f8.1,2i5)') +c & "SSBOND_E_BREAK",totT,t_bath,i,j +#endif +#endif + endif + +c-------TESTING CODE + if (checkstop) then + if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') + & "CHECKSTOP",rij,eij,ed + echeck(jcheck)=eij + endif + enddo + if (checkstop) then + write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps + endif + enddo + if (checkstop) then + transgrad=.true. + checkstop=.false. + endif +c-------END TESTING CODE + + do k=1,3 + dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij + dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij + enddo + do k=1,3 + gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k) + enddo + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k) + & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) + & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + gvdwx(k,j)=gvdwx(k,j)+gg(k) + & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) + & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + enddo +cgrad do k=i,j-1 +cgrad do l=1,3 +cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) +cgrad enddo +cgrad enddo + + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l) + enddo + + return + end + +C----------------------------------------------------------------------------- + + double precision function h_base(x,deriv) +c A smooth function going 0->1 in range [0,1] +c It should NOT be called outside range [0,1], it will not work there. + implicit none + +c Input arguments + double precision x + +c Output arguments + double precision deriv + +c Local variables + double precision xsq + + +c Two parabolas put together. First derivative zero at extrema +c$$$ if (x.lt.0.5D0) then +c$$$ h_base=2.0D0*x*x +c$$$ deriv=4.0D0*x +c$$$ else +c$$$ deriv=1.0D0-x +c$$$ h_base=1.0D0-2.0D0*deriv*deriv +c$$$ deriv=4.0D0*deriv +c$$$ endif + +c Third degree polynomial. First derivative zero at extrema + h_base=x*x*(3.0d0-2.0d0*x) + deriv=6.0d0*x*(1.0d0-x) + +c Fifth degree polynomial. First and second derivatives zero at extrema +c$$$ xsq=x*x +c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0) +c$$$ deriv=x-1.0d0 +c$$$ deriv=deriv*deriv +c$$$ deriv=30.0d0*xsq*deriv + + return + end + +c---------------------------------------------------------------------------- + + subroutine dyn_set_nss +c Adjust nss and other relevant variables based on dyn_ssbond_ij +c implicit none + +c Includes + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" +#endif + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' +#ifndef CLUST + integer maxprocs + parameter (maxprocs=128) +#endif +c#ifdef WHAM + integer max_fg_procs + parameter (max_fg_procs=128) +c#endif +cc include 'COMMON.SETUP' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' + include 'COMMON.SETUP' +#endif +#endif + +c Local variables + double precision emin + integer i,j,imin + integer diff,allflag(maxdim),allnss, + & allihpb(maxdim),alljhpb(maxdim), + & newnss,newihpb(maxdim),newjhpb(maxdim) + logical found + integer i_newnss(max_fg_procs),displ(0:max_fg_procs) + integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss + + allnss=0 + do i=1,nres-1 + do j=i+1,nres + if (dyn_ssbond_ij(i,j).lt.1.0d300) then + allnss=allnss+1 + allflag(allnss)=0 + allihpb(allnss)=i + alljhpb(allnss)=j + endif + enddo + enddo + +cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss) + + 1 emin=1.0d300 + do i=1,allnss + if (allflag(i).eq.0 .and. + & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then + emin=dyn_ssbond_ij(allihpb(i),alljhpb(i)) + imin=i + endif + enddo + if (emin.lt.1.0d300) then + allflag(imin)=1 + do i=1,allnss + if (allflag(i).eq.0 .and. + & (allihpb(i).eq.allihpb(imin) .or. + & alljhpb(i).eq.allihpb(imin) .or. + & allihpb(i).eq.alljhpb(imin) .or. + & alljhpb(i).eq.alljhpb(imin))) then + allflag(i)=-1 + endif + enddo + goto 1 + endif + +cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss) + + newnss=0 + do i=1,allnss + if (allflag(i).eq.1) then + newnss=newnss+1 + newihpb(newnss)=allihpb(i) + newjhpb(newnss)=alljhpb(i) + endif + enddo + +#ifdef MPI + if (nfgtasks.gt.1)then + + call MPI_Reduce(newnss,g_newnss,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) + call MPI_Gather(newnss,1,MPI_INTEGER, + & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_newnss(i-1)+displ(i-1) + enddo + call MPI_Gatherv(newihpb,newnss,MPI_INTEGER, + & g_newihpb,i_newnss,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER, + & g_newjhpb,i_newnss,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + if(fg_rank.eq.0) then +c print *,'g_newnss',g_newnss +c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss) +c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss) + newnss=g_newnss + do i=1,newnss + newihpb(i)=g_newihpb(i) + newjhpb(i)=g_newjhpb(i) + enddo + endif + endif +#endif + + diff=newnss-nss + +cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss) + + do i=1,nss + found=.false. + do j=1,newnss + if (idssb(i).eq.newihpb(j) .and. + & jdssb(i).eq.newjhpb(j)) found=.true. + enddo +#ifndef CLUST +#ifndef WHAM + if (.not.found.and.fg_rank.eq.0) + & write(iout,'(a15,f12.2,f8.1,2i5)') + & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i) +#endif +#endif + enddo + + do i=1,newnss + found=.false. + do j=1,nss + if (newihpb(i).eq.idssb(j) .and. + & newjhpb(i).eq.jdssb(j)) found=.true. + enddo +#ifndef CLUST +#ifndef WHAM + if (.not.found.and.fg_rank.eq.0) + & write(iout,'(a15,f12.2,f8.1,2i5)') + & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) +#endif +#endif + enddo + + nss=newnss + do i=1,nss + idssb(i)=newihpb(i) + jdssb(i)=newjhpb(i) + enddo + + return + end + +c---------------------------------------------------------------------------- +#ifdef NIEWIADOM +c#ifdef WHAM + subroutine read_ssHist + implicit none + +c Includes + include 'DIMENSIONS' + include "DIMENSIONS.FREE" + include 'COMMON.FREE' + +c Local variables + integer i,j + character*80 controlcard + + do i=1,dyn_nssHist + call card_concat(controlcard,.true.) + read(controlcard,*) + & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0)) + enddo + + return + end +#endif + +c---------------------------------------------------------------------------- + + +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- + +c$$$c----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax(i_in,j_in) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c Input arguments +c$$$ integer i_in,j_in +c$$$ +c$$$c Local variables +c$$$ integer i,iretcode,nfun_sc +c$$$ logical scfail +c$$$ double precision var(maxvar),e_sc,etot +c$$$ +c$$$ +c$$$ mask_r=.true. +c$$$ do i=nnt,nct +c$$$ mask_side(i)=0 +c$$$ enddo +c$$$ mask_side(i_in)=1 +c$$$ mask_side(j_in)=1 +c$$$ +c$$$c Minimize the two selected side-chains +c$$$ call overlap_sc(scfail) ! Better not fail! +c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc) +c$$$ +c$$$ mask_r=.false. +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c------------------------------------------------------------- +c$$$ +c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun) +c$$$c Minimize side-chains only, starting from geom but without modifying +c$$$c bond lengths. +c$$$c If mask_r is already set, only the selected side-chains are minimized, +c$$$c otherwise all side-chains are minimized keeping the backbone frozen. +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.MINIM' +c$$$ integer icall +c$$$ common /srutu/ icall +c$$$ +c$$$c Output arguments +c$$$ double precision etot_sc +c$$$ integer iretcode,nfun +c$$$ +c$$$c External functions/subroutines +c$$$ external func_sc,grad_sc,fdum +c$$$ +c$$$c Local variables +c$$$ integer liv,lv +c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) +c$$$ integer iv(liv) +c$$$ double precision rdum(1) +c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar) +c$$$ integer idum(1) +c$$$ integer i,nvar_restr +c$$$ +c$$$ +c$$$cmc start_minim=.true. +c$$$ call deflt(2,iv,liv,lv,v) +c$$$* 12 means fresh start, dont call deflt +c$$$ iv(1)=12 +c$$$* max num of fun calls +c$$$ if (maxfun.eq.0) maxfun=500 +c$$$ iv(17)=maxfun +c$$$* max num of iterations +c$$$ if (maxmin.eq.0) maxmin=1000 +c$$$ iv(18)=maxmin +c$$$* controls output +c$$$ iv(19)=1 +c$$$* selects output unit +c$$$ iv(21)=0 +c$$$c iv(21)=iout ! DEBUG +c$$$c iv(21)=8 ! DEBUG +c$$$* 1 means to print out result +c$$$ iv(22)=0 +c$$$c iv(22)=1 ! DEBUG +c$$$* 1 means to print out summary stats +c$$$ iv(23)=0 +c$$$c iv(23)=1 ! DEBUG +c$$$* 1 means to print initial x and d +c$$$ iv(24)=0 +c$$$c iv(24)=1 ! DEBUG +c$$$* min val for v(radfac) default is 0.1 +c$$$ v(24)=0.1D0 +c$$$* max val for v(radfac) default is 4.0 +c$$$ v(25)=2.0D0 +c$$$c v(25)=4.0D0 +c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +c$$$* the sumsl default is 0.1 +c$$$ v(26)=0.1D0 +c$$$* false conv if (act fnctn decrease) .lt. v(34) +c$$$* the sumsl default is 100*machep +c$$$ v(34)=v(34)/100.0D0 +c$$$* absolute convergence +c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4 +c$$$ v(31)=tolf +c$$$* relative convergence +c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1 +c$$$ v(32)=rtolf +c$$$* controls initial step size +c$$$ v(35)=1.0D-1 +c$$$* large vals of d correspond to small components of step +c$$$ do i=1,nphi +c$$$ d(i)=1.0D-1 +c$$$ enddo +c$$$ do i=nphi+1,nvar +c$$$ d(i)=1.0D-1 +c$$$ enddo +c$$$ +c$$$ call geom_to_var(nvar,x) +c$$$ IF (mask_r) THEN +c$$$ do i=1,nres ! Just in case... +c$$$ mask_phi(i)=0 +c$$$ mask_theta(i)=0 +c$$$ enddo +c$$$ call x2xx(x,xx,nvar_restr) +c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc, +c$$$ & iv,liv,lv,v,idum,rdum,fdum) +c$$$ call xx2x(x,xx) +c$$$ ELSE +c$$$c When minimizing ALL side-chains, etotal_sc is a little +c$$$c faster if we don't set mask_r +c$$$ do i=1,nres +c$$$ mask_phi(i)=0 +c$$$ mask_theta(i)=0 +c$$$ mask_side(i)=1 +c$$$ enddo +c$$$ call x2xx(x,xx,nvar_restr) +c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc, +c$$$ & iv,liv,lv,v,idum,rdum,fdum) +c$$$ call xx2x(x,xx) +c$$$ ENDIF +c$$$ call var_to_geom(nvar,x) +c$$$ call chainbuild_sc +c$$$ etot_sc=v(10) +c$$$ iretcode=iv(1) +c$$$ nfun=iv(6) +c$$$ return +c$$$ end +c$$$ +c$$$C-------------------------------------------------------------------------- +c$$$ +c$$$ subroutine chainbuild_sc +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c Local variables +c$$$ integer i +c$$$ +c$$$ +c$$$ do i=nnt,nct +c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then +c$$$ call locate_side_chain(i) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C-------------------------------------------------------------------------- +c$$$ +c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.MINIM' +c$$$ include 'COMMON.IOUNITS' +c$$$ +c$$$c Input arguments +c$$$ integer n +c$$$ double precision x(maxvar) +c$$$ double precision ufparm +c$$$ external ufparm +c$$$ +c$$$c Input/Output arguments +c$$$ integer nf +c$$$ integer uiparm(1) +c$$$ double precision urparm(1) +c$$$ +c$$$c Output arguments +c$$$ double precision f +c$$$ +c$$$c Local variables +c$$$ double precision energia(0:n_ene) +c$$$#ifdef OSF +c$$$c Variables used to intercept NaNs +c$$$ double precision x_sum +c$$$ integer i_NAN +c$$$#endif +c$$$ +c$$$ +c$$$ nfl=nf +c$$$ icg=mod(nf,2)+1 +c$$$ +c$$$#ifdef OSF +c$$$c Intercept NaNs in the coordinates, before calling etotal_sc +c$$$ x_sum=0.D0 +c$$$ do i_NAN=1,n +c$$$ x_sum=x_sum+x(i_NAN) +c$$$ enddo +c$$$c Calculate the energy only if the coordinates are ok +c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then +c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates" +c$$$ f=1.0D+77 +c$$$ nf=0 +c$$$ else +c$$$#endif +c$$$ +c$$$ call var_to_geom_restr(n,x) +c$$$ call zerograd +c$$$ call chainbuild_sc +c$$$ call etotal_sc(energia(0)) +c$$$ f=energia(0) +c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0 +c$$$ +c$$$#ifdef OSF +c$$$ endif +c$$$#endif +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c------------------------------------------------------- +c$$$ +c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.MINIM' +c$$$ +c$$$c Input arguments +c$$$ integer n +c$$$ double precision x(maxvar) +c$$$ double precision ufparm +c$$$ external ufparm +c$$$ +c$$$c Input/Output arguments +c$$$ integer nf +c$$$ integer uiparm(1) +c$$$ double precision urparm(1) +c$$$ +c$$$c Output arguments +c$$$ double precision g(maxvar) +c$$$ +c$$$c Local variables +c$$$ double precision f,gphii,gthetai,galphai,gomegai +c$$$ integer ig,ind,i,j,k,igall,ij +c$$$ +c$$$ +c$$$ icg=mod(nf,2)+1 +c$$$ if (nf-nfl+1) 20,30,40 +c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm) +c$$$c write (iout,*) 'grad 20' +c$$$ if (nf.eq.0) return +c$$$ goto 40 +c$$$ 30 call var_to_geom_restr(n,x) +c$$$ call chainbuild_sc +c$$$C +c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +c$$$C +c$$$ 40 call cartder +c$$$C +c$$$C Convert the Cartesian gradient into internal-coordinate gradient. +c$$$C +c$$$ +c$$$ ig=0 +c$$$ ind=nres-2 +c$$$ do i=2,nres-2 +c$$$ IF (mask_phi(i+2).eq.1) THEN +c$$$ gphii=0.0D0 +c$$$ do j=i+1,nres-1 +c$$$ ind=ind+1 +c$$$ do k=1,3 +c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) +c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) +c$$$ enddo +c$$$ enddo +c$$$ ig=ig+1 +c$$$ g(ig)=gphii +c$$$ ELSE +c$$$ ind=ind+nres-1-i +c$$$ ENDIF +c$$$ enddo +c$$$ +c$$$ +c$$$ ind=0 +c$$$ do i=1,nres-2 +c$$$ IF (mask_theta(i+2).eq.1) THEN +c$$$ ig=ig+1 +c$$$ gthetai=0.0D0 +c$$$ do j=i+1,nres-1 +c$$$ ind=ind+1 +c$$$ do k=1,3 +c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) +c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) +c$$$ enddo +c$$$ enddo +c$$$ g(ig)=gthetai +c$$$ ELSE +c$$$ ind=ind+nres-1-i +c$$$ ENDIF +c$$$ enddo +c$$$ +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ ig=ig+1 +c$$$ galphai=0.0D0 +c$$$ do k=1,3 +c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg) +c$$$ enddo +c$$$ g(ig)=galphai +c$$$ ENDIF +c$$$ endif +c$$$ enddo +c$$$ +c$$$ +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ ig=ig+1 +c$$$ gomegai=0.0D0 +c$$$ do k=1,3 +c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) +c$$$ enddo +c$$$ g(ig)=gomegai +c$$$ ENDIF +c$$$ endif +c$$$ enddo +c$$$ +c$$$C +c$$$C Add the components corresponding to local energy terms. +c$$$C +c$$$ +c$$$ ig=0 +c$$$ igall=0 +c$$$ do i=4,nres +c$$$ igall=igall+1 +c$$$ if (mask_phi(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ do i=3,nres +c$$$ igall=igall+1 +c$$$ if (mask_theta(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ do ij=1,2 +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ igall=igall+1 +c$$$ if (mask_side(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ endif +c$$$ enddo +c$$$ enddo +c$$$ +c$$$cd do i=1,ig +c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +c$$$cd enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine etotal_sc(energy_sc) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.FFIELD' +c$$$ +c$$$c Output arguments +c$$$ double precision energy_sc(0:n_ene) +c$$$ +c$$$c Local variables +c$$$ double precision evdw,escloc +c$$$ integer i,j +c$$$ +c$$$ +c$$$ do i=1,n_ene +c$$$ energy_sc(i)=0.0D0 +c$$$ enddo +c$$$ +c$$$ if (mask_r) then +c$$$ call egb_sc(evdw) +c$$$ call esc_sc(escloc) +c$$$ else +c$$$ call egb(evdw) +c$$$ call esc(escloc) +c$$$ endif +c$$$ +c$$$ if (evdw.eq.1.0D20) then +c$$$ energy_sc(0)=evdw +c$$$ else +c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc +c$$$ endif +c$$$ energy_sc(1)=evdw +c$$$ energy_sc(12)=escloc +c$$$ +c$$$C +c$$$C Sum up the components of the Cartesian gradient. +c$$$C +c$$$ do i=1,nct +c$$$ do j=1,3 +c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i) +c$$$ enddo +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine egb_sc(evdw) +c$$$C +c$$$C This subroutine calculates the interaction energy of nonbonded side chains +c$$$C assuming the Gay-Berne potential of interaction. +c$$$C +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.CALC' +c$$$ include 'COMMON.CONTROL' +c$$$ logical lprn +c$$$ evdw=0.0D0 +c$$$ energy_dec=.false. +c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +c$$$ evdw=0.0D0 +c$$$ lprn=.false. +c$$$c if (icall.eq.0) lprn=.false. +c$$$ ind=0 +c$$$ do i=iatsc_s,iatsc_e +c$$$ itypi=itype(i) +c$$$ itypi1=itype(i+1) +c$$$ xi=c(1,nres+i) +c$$$ yi=c(2,nres+i) +c$$$ zi=c(3,nres+i) +c$$$ dxi=dc_norm(1,nres+i) +c$$$ dyi=dc_norm(2,nres+i) +c$$$ dzi=dc_norm(3,nres+i) +c$$$c dsci_inv=dsc_inv(itypi) +c$$$ dsci_inv=vbld_inv(i+nres) +c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) +c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi +c$$$C +c$$$C Calculate SC interaction energy. +c$$$C +c$$$ do iint=1,nint_gr(i) +c$$$ do j=istart(i,iint),iend(i,iint) +c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN +c$$$ ind=ind+1 +c$$$ itypj=itype(j) +c$$$c dscj_inv=dsc_inv(itypj) +c$$$ dscj_inv=vbld_inv(j+nres) +c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, +c$$$c & 1.0d0/vbld(j+nres) +c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c$$$ sig0ij=sigma(itypi,itypj) +c$$$ chi1=chi(itypi,itypj) +c$$$ chi2=chi(itypj,itypi) +c$$$ chi12=chi1*chi2 +c$$$ chip1=chip(itypi) +c$$$ chip2=chip(itypj) +c$$$ chip12=chip1*chip2 +c$$$ alf1=alp(itypi) +c$$$ alf2=alp(itypj) +c$$$ alf12=0.5D0*(alf1+alf2) +c$$$C For diagnostics only!!! +c$$$c chi1=0.0D0 +c$$$c chi2=0.0D0 +c$$$c chi12=0.0D0 +c$$$c chip1=0.0D0 +c$$$c chip2=0.0D0 +c$$$c chip12=0.0D0 +c$$$c alf1=0.0D0 +c$$$c alf2=0.0D0 +c$$$c alf12=0.0D0 +c$$$ xj=c(1,nres+j)-xi +c$$$ yj=c(2,nres+j)-yi +c$$$ zj=c(3,nres+j)-zi +c$$$ dxj=dc_norm(1,nres+j) +c$$$ dyj=dc_norm(2,nres+j) +c$$$ dzj=dc_norm(3,nres+j) +c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi +c$$$c write (iout,*) "j",j," dc_norm", +c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c$$$ rij=dsqrt(rrij) +c$$$C Calculate angle-dependent terms of energy and contributions to their +c$$$C derivatives. +c$$$ call sc_angular +c$$$ sigsq=1.0D0/sigsq +c$$$ sig=sig0ij*dsqrt(sigsq) +c$$$ rij_shift=1.0D0/rij-sig+sig0ij +c$$$c for diagnostics; uncomment +c$$$c rij_shift=1.2*sig0ij +c$$$C I hate to put IF's in the loops, but here don't have another choice!!!! +c$$$ if (rij_shift.le.0.0D0) then +c$$$ evdw=1.0D20 +c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$cd & restyp(itypi),i,restyp(itypj),j, +c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) +c$$$ return +c$$$ endif +c$$$ sigder=-sig*sigsq +c$$$c--------------------------------------------------------------- +c$$$ rij_shift=1.0D0/rij_shift +c$$$ fac=rij_shift**expon +c$$$ e1=fac*fac*aa(itypi,itypj) +c$$$ e2=fac*bb(itypi,itypj) +c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2) +c$$$ eps2der=evdwij*eps3rt +c$$$ eps3der=evdwij*eps2rt +c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, +c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +c$$$ evdwij=evdwij*eps2rt*eps3rt +c$$$ evdw=evdw+evdwij +c$$$ if (lprn) then +c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$ & restyp(itypi),i,restyp(itypj),j, +c$$$ & epsi,sigm,chi1,chi2,chip1,chip2, +c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +c$$$ & evdwij +c$$$ endif +c$$$ +c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)') +c$$$ & 'evdw',i,j,evdwij +c$$$ +c$$$C Calculate gradient components. +c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2 +c$$$ fac=-expon*(e1+evdwij)*rij_shift +c$$$ sigder=fac*sigder +c$$$ fac=rij*fac +c$$$c fac=0.0d0 +c$$$C Calculate the radial part of the gradient +c$$$ gg(1)=xj*fac +c$$$ gg(2)=yj*fac +c$$$ gg(3)=zj*fac +c$$$C Calculate angular part of the gradient. +c$$$ call sc_grad +c$$$ ENDIF +c$$$ enddo ! j +c$$$ enddo ! iint +c$$$ enddo ! i +c$$$ energy_dec=.false. +c$$$ return +c$$$ end +c$$$ +c$$$c----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine esc_sc(escloc) +c$$$C Calculate the local energy of a side chain and its derivatives in the +c$$$C corresponding virtual-bond valence angles THETA and the spherical angles +c$$$C ALPHA and OMEGA. +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.FFIELD' +c$$$ include 'COMMON.CONTROL' +c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3), +c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3) +c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit +c$$$ delta=0.02d0*pi +c$$$ escloc=0.0D0 +c$$$c write (iout,'(a)') 'ESC' +c$$$ do i=loc_start,loc_end +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ it=itype(i) +c$$$ if (it.eq.10) goto 1 +c$$$ nlobit=nlob(it) +c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit +c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad +c$$$ theti=theta(i+1)-pipol +c$$$ x(1)=dtan(theti) +c$$$ x(2)=alph(i) +c$$$ x(3)=omeg(i) +c$$$ +c$$$ if (x(2).gt.pi-delta) then +c$$$ xtemp(1)=x(1) +c$$$ xtemp(2)=pi-delta +c$$$ xtemp(3)=x(3) +c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) +c$$$ xtemp(2)=pi +c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.) +c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2), +c$$$ & escloci,dersc(2)) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), +c$$$ & ddersc0(1),dersc(1)) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3), +c$$$ & ddersc0(3),dersc(3)) +c$$$ xtemp(2)=pi-delta +c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) +c$$$ xtemp(2)=pi +c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) +c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1, +c$$$ & dersc0(2),esclocbi,dersc02) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), +c$$$ & dersc12,dersc01) +c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd) +c$$$ dersc0(1)=dersc01 +c$$$ dersc0(2)=dersc02 +c$$$ dersc0(3)=0.0d0 +c$$$ do k=1,3 +c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) +c$$$ enddo +c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c$$$c & esclocbi,ss,ssd +c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi +c$$$c escloci=esclocbi +c$$$c write (iout,*) escloci +c$$$ else if (x(2).lt.delta) then +c$$$ xtemp(1)=x(1) +c$$$ xtemp(2)=delta +c$$$ xtemp(3)=x(3) +c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) +c$$$ xtemp(2)=0.0d0 +c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.) +c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2), +c$$$ & escloci,dersc(2)) +c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), +c$$$ & ddersc0(1),dersc(1)) +c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3), +c$$$ & ddersc0(3),dersc(3)) +c$$$ xtemp(2)=delta +c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) +c$$$ xtemp(2)=0.0d0 +c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) +c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1, +c$$$ & dersc0(2),esclocbi,dersc02) +c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), +c$$$ & dersc12,dersc01) +c$$$ dersc0(1)=dersc01 +c$$$ dersc0(2)=dersc02 +c$$$ dersc0(3)=0.0d0 +c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd) +c$$$ do k=1,3 +c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) +c$$$ enddo +c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c$$$c & esclocbi,ss,ssd +c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi +c$$$c write (iout,*) escloci +c$$$ else +c$$$ call enesc(x,escloci,dersc,ddummy,.false.) +c$$$ endif +c$$$ +c$$$ escloc=escloc+escloci +c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)') +c$$$ & 'escloc',i,escloci +c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc +c$$$ +c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ +c$$$ & wscloc*dersc(1) +c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2) +c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3) +c$$$ 1 continue +c$$$ ENDIF +c$$$ enddo +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine egb_ij(i_sc,j_sc,evdw) +c$$$C +c$$$C This subroutine calculates the interaction energy of nonbonded side chains +c$$$C assuming the Gay-Berne potential of interaction. +c$$$C +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.CALC' +c$$$ include 'COMMON.CONTROL' +c$$$ logical lprn +c$$$ evdw=0.0D0 +c$$$ energy_dec=.false. +c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +c$$$ evdw=0.0D0 +c$$$ lprn=.false. +c$$$ ind=0 +c$$$c$$$ do i=iatsc_s,iatsc_e +c$$$ i=i_sc +c$$$ itypi=itype(i) +c$$$ itypi1=itype(i+1) +c$$$ xi=c(1,nres+i) +c$$$ yi=c(2,nres+i) +c$$$ zi=c(3,nres+i) +c$$$ dxi=dc_norm(1,nres+i) +c$$$ dyi=dc_norm(2,nres+i) +c$$$ dzi=dc_norm(3,nres+i) +c$$$c dsci_inv=dsc_inv(itypi) +c$$$ dsci_inv=vbld_inv(i+nres) +c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) +c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi +c$$$C +c$$$C Calculate SC interaction energy. +c$$$C +c$$$c$$$ do iint=1,nint_gr(i) +c$$$c$$$ do j=istart(i,iint),iend(i,iint) +c$$$ j=j_sc +c$$$ ind=ind+1 +c$$$ itypj=itype(j) +c$$$c dscj_inv=dsc_inv(itypj) +c$$$ dscj_inv=vbld_inv(j+nres) +c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, +c$$$c & 1.0d0/vbld(j+nres) +c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c$$$ sig0ij=sigma(itypi,itypj) +c$$$ chi1=chi(itypi,itypj) +c$$$ chi2=chi(itypj,itypi) +c$$$ chi12=chi1*chi2 +c$$$ chip1=chip(itypi) +c$$$ chip2=chip(itypj) +c$$$ chip12=chip1*chip2 +c$$$ alf1=alp(itypi) +c$$$ alf2=alp(itypj) +c$$$ alf12=0.5D0*(alf1+alf2) +c$$$C For diagnostics only!!! +c$$$c chi1=0.0D0 +c$$$c chi2=0.0D0 +c$$$c chi12=0.0D0 +c$$$c chip1=0.0D0 +c$$$c chip2=0.0D0 +c$$$c chip12=0.0D0 +c$$$c alf1=0.0D0 +c$$$c alf2=0.0D0 +c$$$c alf12=0.0D0 +c$$$ xj=c(1,nres+j)-xi +c$$$ yj=c(2,nres+j)-yi +c$$$ zj=c(3,nres+j)-zi +c$$$ dxj=dc_norm(1,nres+j) +c$$$ dyj=dc_norm(2,nres+j) +c$$$ dzj=dc_norm(3,nres+j) +c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi +c$$$c write (iout,*) "j",j," dc_norm", +c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c$$$ rij=dsqrt(rrij) +c$$$C Calculate angle-dependent terms of energy and contributions to their +c$$$C derivatives. +c$$$ call sc_angular +c$$$ sigsq=1.0D0/sigsq +c$$$ sig=sig0ij*dsqrt(sigsq) +c$$$ rij_shift=1.0D0/rij-sig+sig0ij +c$$$c for diagnostics; uncomment +c$$$c rij_shift=1.2*sig0ij +c$$$C I hate to put IF's in the loops, but here don't have another choice!!!! +c$$$ if (rij_shift.le.0.0D0) then +c$$$ evdw=1.0D20 +c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$cd & restyp(itypi),i,restyp(itypj),j, +c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) +c$$$ return +c$$$ endif +c$$$ sigder=-sig*sigsq +c$$$c--------------------------------------------------------------- +c$$$ rij_shift=1.0D0/rij_shift +c$$$ fac=rij_shift**expon +c$$$ e1=fac*fac*aa(itypi,itypj) +c$$$ e2=fac*bb(itypi,itypj) +c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2) +c$$$ eps2der=evdwij*eps3rt +c$$$ eps3der=evdwij*eps2rt +c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, +c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +c$$$ evdwij=evdwij*eps2rt*eps3rt +c$$$ evdw=evdw+evdwij +c$$$ if (lprn) then +c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$ & restyp(itypi),i,restyp(itypj),j, +c$$$ & epsi,sigm,chi1,chi2,chip1,chip2, +c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +c$$$ & evdwij +c$$$ endif +c$$$ +c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)') +c$$$ & 'evdw',i,j,evdwij +c$$$ +c$$$C Calculate gradient components. +c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2 +c$$$ fac=-expon*(e1+evdwij)*rij_shift +c$$$ sigder=fac*sigder +c$$$ fac=rij*fac +c$$$c fac=0.0d0 +c$$$C Calculate the radial part of the gradient +c$$$ gg(1)=xj*fac +c$$$ gg(2)=yj*fac +c$$$ gg(3)=zj*fac +c$$$C Calculate angular part of the gradient. +c$$$ call sc_grad +c$$$c$$$ enddo ! j +c$$$c$$$ enddo ! iint +c$$$c$$$ enddo ! i +c$$$ energy_dec=.false. +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine perturb_side_chain(i,angle) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.IOUNITS' +c$$$ +c$$$c External functions +c$$$ external ran_number +c$$$ double precision ran_number +c$$$ +c$$$c Input arguments +c$$$ integer i +c$$$ double precision angle ! In degrees +c$$$ +c$$$c Local variables +c$$$ integer i_sc +c$$$ double precision rad_ang,rand_v(3),length,cost,sint +c$$$ +c$$$ +c$$$ i_sc=i+nres +c$$$ rad_ang=angle*deg2rad +c$$$ +c$$$ length=0.0 +c$$$ do while (length.lt.0.01) +c$$$ rand_v(1)=ran_number(0.01D0,1.0D0) +c$$$ rand_v(2)=ran_number(0.01D0,1.0D0) +c$$$ rand_v(3)=ran_number(0.01D0,1.0D0) +c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+ +c$$$ + rand_v(3)*rand_v(3) +c$$$ length=sqrt(length) +c$$$ rand_v(1)=rand_v(1)/length +c$$$ rand_v(2)=rand_v(2)/length +c$$$ rand_v(3)=rand_v(3)/length +c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+ +c$$$ + rand_v(3)*dc_norm(3,i_sc) +c$$$ length=1.0D0-cost*cost +c$$$ if (length.lt.0.0D0) length=0.0D0 +c$$$ length=sqrt(length) +c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc) +c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc) +c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc) +c$$$ enddo +c$$$ rand_v(1)=rand_v(1)/length +c$$$ rand_v(2)=rand_v(2)/length +c$$$ rand_v(3)=rand_v(3)/length +c$$$ +c$$$ cost=dcos(rad_ang) +c$$$ sint=dsin(rad_ang) +c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint) +c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint) +c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint) +c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc) +c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc) +c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc) +c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc) +c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc) +c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc) +c$$$ +c$$$ call chainbuild_cart +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c---------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax3(i_in,j_in) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c External functions +c$$$ external ran_number +c$$$ double precision ran_number +c$$$ +c$$$c Input arguments +c$$$ integer i_in,j_in +c$$$ +c$$$c Local variables +c$$$ double precision energy_sc(0:n_ene),etot +c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3) +c$$$ double precision ang_pert,rand_fact,exp_fact,beta +c$$$ integer n,i_pert,i +c$$$ logical notdone +c$$$ +c$$$ +c$$$ beta=1.0D0 +c$$$ +c$$$ mask_r=.true. +c$$$ do i=nnt,nct +c$$$ mask_side(i)=0 +c$$$ enddo +c$$$ mask_side(i_in)=1 +c$$$ mask_side(j_in)=1 +c$$$ +c$$$ call etotal_sc(energy_sc) +c$$$ etot=energy_sc(0) +c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ +c$$$ notdone=.true. +c$$$ n=0 +c$$$ do while (notdone) +c$$$ if (mod(n,2).eq.0) then +c$$$ i_pert=i_in +c$$$ else +c$$$ i_pert=j_in +c$$$ endif +c$$$ n=n+1 +c$$$ +c$$$ do i=1,3 +c$$$ org_dc(i)=dc(i,i_pert+nres) +c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres) +c$$$ org_c(i)=c(i,i_pert+nres) +c$$$ enddo +c$$$ ang_pert=ran_number(0.0D0,3.0D0) +c$$$ call perturb_side_chain(i_pert,ang_pert) +c$$$ call etotal_sc(energy_sc) +c$$$ exp_fact=exp(beta*(etot-energy_sc(0))) +c$$$ rand_fact=ran_number(0.0D0,1.0D0) +c$$$ if (rand_fact.lt.exp_fact) then +c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ etot=energy_sc(0) +c$$$ else +c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ do i=1,3 +c$$$ dc(i,i_pert+nres)=org_dc(i) +c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i) +c$$$ c(i,i_pert+nres)=org_c(i) +c$$$ enddo +c$$$ endif +c$$$ +c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false. +c$$$ enddo +c$$$ +c$$$ mask_r=.false. +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c---------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in) +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ integer liv,lv +c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2)) +c$$$********************************************************************* +c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for * +c$$$* the calling subprogram. * +c$$$* when d(i)=1.0, then v(35) is the length of the initial step, * +c$$$* calculated in the usual pythagorean way. * +c$$$* absolute convergence occurs when the function is within v(31) of * +c$$$* zero. unless you know the minimum value in advance, abs convg * +c$$$* is probably not useful. * +c$$$* relative convergence is when the model predicts that the function * +c$$$* will decrease by less than v(32)*abs(fun). * +c$$$********************************************************************* +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.MINIM' +c$$$ include 'COMMON.CHAIN' +c$$$ +c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist +c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar), +c$$$ + orig_ss_dist(maxres2,maxres2) +c$$$ +c$$$ double precision etot +c$$$ integer iretcode,nfun,i_in,j_in +c$$$ +c$$$ external dist +c$$$ double precision dist +c$$$ external ss_func,fdum +c$$$ double precision ss_func,fdum +c$$$ +c$$$ integer iv(liv),uiparm(2) +c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum +c$$$ integer i,j,k +c$$$ +c$$$ +c$$$ call deflt(2,iv,liv,lv,v) +c$$$* 12 means fresh start, dont call deflt +c$$$ iv(1)=12 +c$$$* max num of fun calls +c$$$ if (maxfun.eq.0) maxfun=500 +c$$$ iv(17)=maxfun +c$$$* max num of iterations +c$$$ if (maxmin.eq.0) maxmin=1000 +c$$$ iv(18)=maxmin +c$$$* controls output +c$$$ iv(19)=2 +c$$$* selects output unit +c$$$c iv(21)=iout +c$$$ iv(21)=0 +c$$$* 1 means to print out result +c$$$ iv(22)=0 +c$$$* 1 means to print out summary stats +c$$$ iv(23)=0 +c$$$* 1 means to print initial x and d +c$$$ iv(24)=0 +c$$$* min val for v(radfac) default is 0.1 +c$$$ v(24)=0.1D0 +c$$$* max val for v(radfac) default is 4.0 +c$$$ v(25)=2.0D0 +c$$$c v(25)=4.0D0 +c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +c$$$* the sumsl default is 0.1 +c$$$ v(26)=0.1D0 +c$$$* false conv if (act fnctn decrease) .lt. v(34) +c$$$* the sumsl default is 100*machep +c$$$ v(34)=v(34)/100.0D0 +c$$$* absolute convergence +c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4 +c$$$ v(31)=tolf +c$$$ v(31)=1.0D-1 +c$$$* relative convergence +c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4 +c$$$ v(32)=rtolf +c$$$ v(32)=1.0D-1 +c$$$* controls initial step size +c$$$ v(35)=1.0D-1 +c$$$* large vals of d correspond to small components of step +c$$$ do i=1,6*nres +c$$$ d(i)=1.0D0 +c$$$ enddo +c$$$ +c$$$ do i=0,2*nres +c$$$ do j=1,3 +c$$$ orig_ss_dc(j,i)=dc(j,i) +c$$$ enddo +c$$$ enddo +c$$$ call geom_to_var(nvar,orig_ss_var) +c$$$ +c$$$ do i=1,nres +c$$$ do j=i,nres +c$$$ orig_ss_dist(j,i)=dist(j,i) +c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i) +c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres) +c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres) +c$$$ enddo +c$$$ enddo +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ x(k)=dc(j,i) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ x(k)=dc(j,i+nres) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ +c$$$ uiparm(1)=i_in +c$$$ uiparm(2)=j_in +c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum) +c$$$ etot=v(10) +c$$$ iretcode=iv(1) +c$$$ nfun=iv(6)+iv(30) +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i)=x(k) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i+nres)=x(k) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ call chainbuild_cart +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.SBRIDGE' +c$$$ +c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist +c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar), +c$$$ + orig_ss_dist(maxres2,maxres2) +c$$$ +c$$$ integer n +c$$$ double precision x(maxres6) +c$$$ integer nf +c$$$ double precision f +c$$$ integer uiparm(2) +c$$$ real*8 urparm(1) +c$$$ external ufparm +c$$$ double precision ufparm +c$$$ +c$$$ external dist +c$$$ double precision dist +c$$$ +c$$$ integer i,j,k,ss_i,ss_j +c$$$ double precision tempf,var(maxvar) +c$$$ +c$$$ +c$$$ ss_i=uiparm(1) +c$$$ ss_j=uiparm(2) +c$$$ f=0.0D0 +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i)=x(k) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i+nres)=x(k) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ call chainbuild_cart +c$$$ +c$$$ call geom_to_var(nvar,var) +c$$$ +c$$$c Constraints on all angles +c$$$ do i=1,nvar +c$$$ tempf=var(i)-orig_ss_var(i) +c$$$ f=f+tempf*tempf +c$$$ enddo +c$$$ +c$$$c Constraints on all distances +c$$$ do i=1,nres-1 +c$$$ if (i.gt.1) then +c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i) +c$$$ f=f+tempf*tempf +c$$$ endif +c$$$ do j=i+1,nres +c$$$ tempf=dist(j,i)-orig_ss_dist(j,i) +c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf +c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ enddo +c$$$ enddo +c$$$ +c$$$c Constraints for the relevant CYS-CYS +c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0 +c$$$ f=f+tempf*tempf +c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF +c$$$ +c$$$c$$$ if (nf.ne.nfl) then +c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf, +c$$$c$$$ + f,dist(5+nres,14+nres) +c$$$c$$$ endif +c$$$ +c$$$ nfl=nf +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- diff --git a/source/wham/src/COMMON.ALLPARM b/source/wham/src/COMMON.ALLPARM index 896b5a2..0dbe565 100644 --- a/source/wham/src/COMMON.ALLPARM +++ b/source/wham/src/COMMON.ALLPARM @@ -47,6 +47,7 @@ & sigma_all(ntyp,ntyp,max_parm),r0_all(ntyp,ntyp,max_parm), & chi_all(ntyp,ntyp,max_parm),chip_all(ntyp,max_parm), & alp_all(ntyp,max_parm),ebr_all(max_parm),d0cm_all(max_parm), + & ss_depth_all(max_parm),ht_all(max_parm), & akcm_all(max_parm),akth_all(max_parm),akct_all(max_parm), & v1ss_all(max_parm),v2ss_all(max_parm),v3ss_all(max_parm), & v1sccor_all(maxterm_sccor,3,ntyp,ntyp,max_parm), @@ -71,6 +72,7 @@ & ctilde_all,dtilde_all,b1tilde_all,app_all,bpp_all,ael6_all, & ael3_all,aad_all,bad_all,aa_all,bb_all,augm_all, & eps_all,sigma_all,r0_all,chi_all,chip_all,alp_all,ebr_all, + & ss_depth_all,ht_all, & d0cm_all,akcm_all,akth_all,akct_all,v1ss_all,v2ss_all,v3ss_all, & v1sccor_all,v2sccor_all,nbondterm_all, & nlob_all,nlor_all,nterm_all,ntermd1_all,ntermd2_all, diff --git a/source/wham/src/bxread.F b/source/wham/src/bxread.F index 32c935d..bebf420 100644 --- a/source/wham/src/bxread.F +++ b/source/wham/src/bxread.F @@ -48,6 +48,8 @@ & eini,efree,rmsdev,(prop(j),j=1,nQ),iscor ii=ii+1 kk=kk+1 + write(iout,*) 'BXWEJ',eini,l + flush(iout) if (mod(kk,isampl(iparm)).eq.0) then jj=jj+1 write(ientout,rec=jj) diff --git a/source/wham/src/cxread.F b/source/wham/src/cxread.F index 0e6c52e..a662f7a 100644 --- a/source/wham/src/cxread.F +++ b/source/wham/src/cxread.F @@ -71,8 +71,10 @@ c print *,"bumbum" #endif do j=1,nss if (dyn_ss) then - call xdrfint_(ixdrf, idssb(j)+nres, iret) - call xdrfint_(ixdrf, jdssb(j)+nres, iret) + call xdrfint_(ixdrf, idssb(j), iret) + call xdrfint_(ixdrf, jdssb(j), iret) + idssb(j)=idssb(j)-nres + jdssb(j)=jdssb(j)-nres else call xdrfint_(ixdrf, ihpb(j), iret) call xdrfint_(ixdrf, jhpb(j), iret) @@ -102,8 +104,11 @@ c print *,"bumbum" #endif do j=1,nss if (dyn_ss) then - call xdrfint(ixdrf, idssb(j)+nres, iret) - call xdrfint(ixdrf, jdssb(j)+nres, iret) + call xdrfint(ixdrf, idssb(j), iret) + call xdrfint(ixdrf, jdssb(j), iret) +cc idssb(j)=idssb(j)-nres +cc jdssb(j)=jdssb(j)-nres +cc write(iout,*) idssb(j),jdssb(j) else call xdrfint(ixdrf, ihpb(j), iret) call xdrfint(ixdrf, jhpb(j), iret) diff --git a/source/wham/src/enecalc1.F b/source/wham/src/enecalc1.F index 34a5eec..01e5684 100644 --- a/source/wham/src/enecalc1.F +++ b/source/wham/src/enecalc1.F @@ -75,6 +75,7 @@ & ((csingle(l,k+nres),l=1,3),k=nnt,nct), & nss,(ihpb(k),jhpb(k),k=1,nss), & eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar +cc write(iout,*), 'NAWEJ',i,eini if (indpdb.gt.0) then do k=1,nres do l=1,3 @@ -217,6 +218,8 @@ c call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev) do k=1,21 enetb(k,iii+1,iparm)=energia(k) enddo +c write(iout,*) "iCHUJ TU STRZELI",i,iii,entfac(i) +c call enerprint(energia(0),fT) #ifdef DEBUG write (iout,'(2i5,f10.1,3e15.5)') i,iii, & 1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree @@ -440,12 +443,22 @@ c------------------------------------------------------------------------------ #else do i=1,ntot(islice) #endif +cc if (dyn_ss) then +cc read(ientout,rec=i,err=101) +cc & ((csingle(l,k),l=1,3),k=1,nres), +cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct), +cc & nss,(idssb(k),jdssb(k),k=1,nss), +cc & eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm +cc idssb(k)=idssb(k)-nres +cc jdssb(k)=jdssb(k)-nres +cc else 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), & nss,(ihpb(k),jhpb(k),k=1,nss), & eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm -c write (iout,*) iR,ib,iparm,eini,efree +cc endif +cc write (iout,*) 'CC', iR,ib,iparm,eini,efree do j=1,2*nres do k=1,3 c(k,j)=csingle(k,j) @@ -455,14 +468,24 @@ c write (iout,*) iR,ib,iparm,eini,efree iscore=0 if (indpdb.gt.0) then call conf_compar(i,.false.,.true.) - endif + endif +c if (dyn_ss) then if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) & ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k+nres),l=1,3),k=nnt,nct), & nss,(ihpb(k),jhpb(k),k=1,nss), c & potE(i,iparm),-entfac(i),rms_nat,iscore & potE(i,nparmset),-entfac(i),rms_nat,iscore -c write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i) +c else + if (bxfile .or.cxfile .or. ensembles.gt.0) write + & (ientin,rec=i) + & ((csingle(l,k),l=1,3),k=1,nres), + & ((csingle(l,k+nres),l=1,3),k=nnt,nct), + & nss,(ihpb(k),jhpb(k),k=1,nss), +c & potE(i,iparm),-entfac(i),rms_nat,iscore + & potE(i,nparmset),-entfac(i),rms_nat,iscore +c endif + write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i) #ifndef MPI if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset), & -entfac(i),rms_nat,iscore) @@ -541,17 +564,37 @@ c write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j) c call flush(iout) do i=indstart(j),indend(j) iii = iii+1 +cc if (dyn_ss) then +cc read(ientin,rec=iii,err=101) +cc & ((csingle(l,k),l=1,3),k=1,nres), +cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct), +cc & nss,(idssb(k),jdssb(k),k=1,nss), +cc & eini,efree,rmsdev,iscor +cc idssb(k)=idssb(k)-nres +cc jdssb(k)=jdssb(k)-nres +cc else read(ientin,rec=iii,err=101) & ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k+nres),l=1,3),k=nnt,nct), & nss,(ihpb(k),jhpb(k),k=1,nss), & eini,efree,rmsdev,iscor +cc endif if (bxfile .or. ensembles.gt.0) then - write (ientout,rec=i) +cc if (dyn_ss) then +cc write (ientout,rec=i) +cc & ((csingle(l,k),l=1,3),k=1,nres), +cc & ((csingle(l,k+nres),l=1,3),k=nnt,nct), +cc & nss,(idssb(k)+nres,jdssb(k)+nres,k=1,nss), +cc & eini,efree,rmsdev,iscor +cc else + write (ientout,rec=i) & ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k+nres),l=1,3),k=nnt,nct), & nss,(ihpb(k),jhpb(k),k=1,nss), & eini,efree,rmsdev,iscor +cc write(iout,*) "W poszukiwaniu zlotych galotow" +cc write(iout,*) "efree=",efree,iii +cc endif endif if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor) #ifdef DEBUG @@ -645,6 +688,7 @@ c call flush(iout) c write (iout,*) "itmp",itmp c call flush(iout) +c write (iout,*) "CNZ",eini,dyn_ss #if (defined(AIX) && !defined(JUBL)) call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret) @@ -652,8 +696,13 @@ c write (iout,*) "xdrf3dfcoord" c call flush(iout) call xdrfint_(ixdrf, nss, iret) do j=1,nss +cc if (dyn_ss) then +cc call xdrfint_(ixdrf, idssb(j)+nres, iret) +cc call xdrfint_(ixdrf, jdssb(j)+nres, iret) +cc else call xdrfint_(ixdrf, ihpb(j), iret) call xdrfint_(ixdrf, jhpb(j), iret) +cc endif enddo call xdrffloat_(ixdrf,real(eini),iret) call xdrffloat_(ixdrf,real(efree),iret) @@ -664,11 +713,18 @@ c call flush(iout) call xdrfint(ixdrf, nss, iret) do j=1,nss +cc if (dyn_ss) then +cc call xdrfint(ixdrf, idssb(j), iret) +cc call xdrfint(ixdrf, jdssb(j), iret) +cc idssb(j)=idssb(j)-nres +cc jdssb(j)=jdssb(j)-nres +cc else call xdrfint(ixdrf, ihpb(j), iret) call xdrfint(ixdrf, jhpb(j), iret) +cc endif enddo call xdrffloat(ixdrf,real(eini),iret) - call xdrffloat(ixdrf,real(efree),iret) + call xdrffloat(ixdrf,real(efree),iret) call xdrffloat(ixdrf,real(rmsdev),iret) call xdrfint(ixdrf,iscor,iret) #endif diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index 9ea237d..9e0b2c5 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -154,7 +154,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t - if (dyn_ss) call dyn_set_nss +c if (dyn_ss) call dyn_set_nss c detecting NaNQ #ifdef ISNAN #ifdef AIX diff --git a/source/wham/src/geomout.F b/source/wham/src/geomout.F index 9506eed..35f3d77 100644 --- a/source/wham/src/geomout.F +++ b/source/wham/src/geomout.F @@ -40,7 +40,11 @@ write (ipdb,30) ica(nct),ica(nct)+1 endif do i=1,nss +c if(dyn_ss) then +c write (ipdb,30) ica(idssb(i))+1,ica(jdssb(i))+1 +c else write (ipdb,30) ica(ihpb(i))+1,ica(jhpb(i))+1 +c endif enddo 10 FORMAT ('ATOM',I7,' CA ',A3,I6,4X,3F8.3) 20 FORMAT ('ATOM',I7,' CB ',A3,I6,4X,3F8.3) diff --git a/source/wham/src/make_ensemble1.F b/source/wham/src/make_ensemble1.F index 5d7b750..e9c0754 100644 --- a/source/wham/src/make_ensemble1.F +++ b/source/wham/src/make_ensemble1.F @@ -166,7 +166,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc @@ -176,7 +176,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr @@ -188,7 +188,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & beta_h(ib,iparm)*etot-entfac(i) potE(i,iparm)=etot #ifdef DEBUG - write (iout,*) i,indstart(me)+i-1,ib, + write (iout,*) 'EEE',i,indstart(me)+i-1,ib, & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm), & -entfac(i),Fdimless(i) #endif @@ -336,7 +336,15 @@ c write (iout,*) "qfree",qfree enddo eini=fdimless(i) call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev) +cc if (temper.eq.300.0d0) then +cc write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree +cc flush(iout) +cc endif enddo +cc if (temper.eq.300.0d0) then +cc write(iout,*) 'Gral i=',i ,eini,enepot(i),efree +cc flush(iout) +cc endif #ifdef MPI endif #endif diff --git a/source/wham/src/parmread.F b/source/wham/src/parmread.F index aae7c01..36730f5 100644 --- a/source/wham/src/parmread.F +++ b/source/wham/src/parmread.F @@ -961,7 +961,7 @@ C C C Define the constants of the disulfide bridge C - ebr=-5.50D0 +c ebr=-5.50D0 c c Old arbitrary potential - commented out. c @@ -972,21 +972,21 @@ c Constants of the disulfide-bond potential determined based on the RHF/6-31G** c energy surface of diethyl disulfide. c A. Liwo and U. Kozlowska, 11/24/03 c - D0CM = 3.78d0 - AKCM = 15.1d0 - AKTH = 11.0d0 - AKCT = 12.0d0 - V1SS =-1.08d0 - V2SS = 7.61d0 - V3SS = 13.7d0 +c D0CM = 3.78d0 +c AKCM = 15.1d0 +c AKTH = 11.0d0 +c AKCT = 12.0d0 +c V1SS =-1.08d0 +c V2SS = 7.61d0 +c V3SS = 13.7d0 - if (lprint) then - write (iout,'(/a)') "Disulfide bridge parameters:" - write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr - write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm - write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct - write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, - & ' v3ss:',v3ss - endif +c if (lprint) then +c write (iout,'(/a)') "Disulfide bridge parameters:" +c write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr +c write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm +c write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct +c write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, +c & ' v3ss:',v3ss +c endif return end diff --git a/source/wham/src/store_parm.F b/source/wham/src/store_parm.F index 6aa33c5..ce0047f 100644 --- a/source/wham/src/store_parm.F +++ b/source/wham/src/store_parm.F @@ -217,6 +217,8 @@ c Store the SCp parameters enddo enddo c Store disulfide-bond parameters + ht_all(iparm)=ht + ss_depth_all(iparm)=ss_depth ebr_all(iparm)=ebr d0cm_all(iparm)=d0cm akcm_all(iparm)=akcm @@ -460,6 +462,8 @@ c Restore the SCp parameters enddo enddo c Restore disulfide-bond parameters + ht=ht_all(iparm) + ss_depth=ss_depth_all(iparm) ebr=ebr_all(iparm) d0cm=d0cm_all(iparm) akcm=akcm_all(iparm) diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F index a6044cd..676d4f4 100644 --- a/source/wham/src/wham_calc1.F +++ b/source/wham/src/wham_calc1.F @@ -316,7 +316,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc @@ -326,7 +326,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr @@ -675,7 +675,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc @@ -685,7 +685,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr @@ -1013,7 +1013,7 @@ c write (iout,*) ib," PotEmin",potEmin etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc @@ -1036,7 +1036,7 @@ c write (iout,*) ib," PotEmin",potEmin etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr -- 1.7.9.5