call readi(controlcard,'CONSTR_DIST',constr_dist,0)
write (iout,*) "with_dihed_constr ",with_dihed_constr,
& " CONSTR_DIST",constr_dist
-<<<<<<< HEAD
-=======
-
- call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
- write (iout,*) "with_homology_constr ",with_dihed_constr,
- & " CONSTR_HOMOLOGY",constr_homology
- print_homology_restraints=
- & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
- print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
- print_homology_models=
- & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
-
-#ifdef AIX
- call flush_(iout)
-#else
->>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit
call flush(iout)
if (min_var) iopt=1
return
if (nend.gt.nct .or. nend.eq.0) nend=nct
write (iout,*) "nstart",nstart," nend",nend
nres0=nres
-<<<<<<< HEAD
-=======
-
-C Juyong:READ init_vars
-C Initialize variables!
-C Juyong:READ read_info
-C READ fragment information!!
-C both routines should be in dfa.F file!!
-
- if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
- & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
-#ifdef DEBUG
- write (iout,*) "Calling init_dfa_vars"
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-#endif
- call init_dfa_vars
-#ifdef DEBUG
- write (iout,*) 'init_dfa_vars finished!'
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-#endif
- call read_dfa_info
-#ifdef DEBUG
- write (iout,*) 'read_dfa_info finished!'
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-#endif
- endif
-
- if (constr_homology.gt.0) then
- call read_constr_homology(print_homology_restraints)
- endif
-
->>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit
c if (pdbref) then
c read(inp,'(a)') pdbfile
c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
call flush(iout)
return
end
-<<<<<<< HEAD
-=======
-
-c====-------------------------------------------------------------------
- subroutine read_constr_homology(lprn)
-
- include 'DIMENSIONS'
-#ifdef MPI
- include 'mpif.h'
-#endif
- include 'COMMON.SETUP'
- include 'COMMON.CONTROL'
- include 'COMMON.CHAIN'
- include 'COMMON.IOUNITS'
- include 'COMMON.GEO'
- include 'COMMON.INTERACT'
- include 'COMMON.HOMRESTR'
-c
-c For new homol impl
-c
- include 'COMMON.VAR'
-c include 'include_unres/COMMON.VAR'
-c
-
-c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
-c & dist_cut
-c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
-c & sigma_odl_temp(maxres,maxres,max_template)
- character*2 kic2
- character*24 model_ki_dist, model_ki_angle
- character*500 controlcard
- integer ki, i, j, k, l
- logical lprn
- logical unres_pdb
-c
-c FP - Nov. 2014 Temporary specifications for new vars
-c
- double precision rescore_tmp,x12,y12,z12
- double precision, dimension (max_template,maxres) :: rescore
- character*24 tpl_k_rescore
-c -----------------------------------------------------------------
-c Reading multiple PDB ref structures and calculation of retraints
-c not using pre-computed ones stored in files model_ki_{dist,angle}
-c FP (Nov., 2014)
-c -----------------------------------------------------------------
-c
-c
-c Alternative: reading from input
-#ifdef DEBUG
- write (iout,*) "BEGIN READ HOMOLOGY INFO"
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-#endif
- call card_concat(controlcard)
- call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
- call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
- call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
- call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
- call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
-
- call readi(controlcard,"HOMOL_NSET",homol_nset,1)
- if (homol_nset.gt.1)then
- call readi(controlcard,"ISET",iset,1)
- call card_concat(controlcard)
- read(controlcard,*) (waga_homology(i),i=1,homol_nset)
- else
- iset=1
- waga_homology(1)=1.0
- endif
-c
-#ifdef DEBUG
- write(iout,*) "read_constr_homology iset",iset
- write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-#endif
-cd write (iout,*) "nnt",nnt," nct",nct
-cd call flush(iout)
-
-
- lim_odl=0
- lim_dih=0
-c
-c New
-c
- lim_theta=0
- lim_xx=0
-c
-c Reading HM global scores (prob not required)
-c
-c open (4,file="HMscore")
-c do k=1,constr_homology
-c read (4,*,end=521) hmscore_tmp
-c hmscore(k)=hmscore_tmp ! Another transformation can be used
-c write(*,*) "Model", k, ":", hmscore(k)
-c enddo
-c521 continue
-
-c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
-
- write (iout,*) "CONSTR_HOMOLOGY",constr_homology
- do k=1,constr_homology
-
- read(inp,'(a)') pdbfile
-c write (iout,*) "k ",k," pdbfile ",pdbfile
-c Next stament causes error upon compilation (?)
-c if(me.eq.king.or. .not. out1file)
-c write (iout,'(2a)') 'PDB data will be read from file ',
-c & pdbfile(:ilen(pdbfile))
- open(ipdbin,file=pdbfile,status='old',err=33)
- goto 34
- 33 write (iout,'(a)') 'Error opening PDB file.'
- stop
- 34 continue
-c print *,'Begin reading pdb data'
-c
-c Files containing res sim or local scores (former containing sigmas)
-c
-
- write(kic2,'(bz,i2.2)') k
-
- tpl_k_rescore="template"//kic2//".sco"
-c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
-c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
-c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
-c tpl_k_sigma_d="template"//kic2//".sigma_d"
-
- unres_pdb=.false.
- call readpdb
- do i=1,2*nres
- do j=1,3
- crefjlee(j,i)=c(j,i)
- enddo
- enddo
-#ifdef DEBUG
- do i=1,nres
- write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
- & (crefjlee(j,i+nres),j=1,3)
- enddo
- write (iout,*) "READ HOMOLOGY INFO"
- write (iout,*) "read_constr_homology x: after reading pdb file"
- write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
- write (iout,*) "waga_dist",waga_dist
- write (iout,*) "waga_angle",waga_angle
- write (iout,*) "waga_theta",waga_theta
- write (iout,*) "waga_d",waga_d
- write (iout,*) "dist_cut",dist_cut
-#endif
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
-
-c
-c Distance restraints
-c
-c ... --> odl(k,ii)
-C Copy the coordinates from reference coordinates (?)
- do i=1,2*nres
- do j=1,3
- c(j,i)=cref(j,i)
-c write (iout,*) "c(",j,i,") =",c(j,i)
- enddo
- enddo
-c
-c From read_dist_constr (commented out 25/11/2014 <-> res sim)
-c
-c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
- open (ientin,file=tpl_k_rescore,status='old')
- do irec=1,maxdim ! loop for reading res sim
- if (irec.eq.1) then
- rescore(k,irec)=0.0d0
- goto 1301
- endif
- read (ientin,*,end=1401) rescore_tmp
-c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
- rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
-c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
- 1301 continue
- enddo
- 1401 continue
- close (ientin)
-c open (ientin,file=tpl_k_sigma_odl,status='old')
-c do irec=1,maxdim ! loop for reading sigma_odl
-c read (ientin,*,end=1401) i, j,
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
-c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
-c enddo
-c 1401 continue
-c close (ientin)
- if (waga_dist.ne.0.0d0) then
- ii=0
- do i = nnt,nct-2 ! right? without parallel.
- do j=i+2,nct ! right?
-c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do j=i+2,nres ! ibid
-c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
-c do j=i+2,nct ! ibid
- ii=ii+1
-c write (iout,*) "k",k
-c write (iout,*) "i",i," j",j," constr_homology",
-c & constr_homology
- ires_homo(ii)=i
- jres_homo(ii)=j
-c
-c Attempt to replace dist(i,j) by its definition in ...
-c
- x12=c(1,i)-c(1,j)
- y12=c(2,i)-c(2,j)
- z12=c(3,i)-c(3,j)
- distal=dsqrt(x12*x12+y12*y12+z12*z12)
- odl(k,ii)=distal
-c
-c odl(k,ii)=dist(i,j)
-c write (iout,*) "dist(",i,j,") =",dist(i,j)
-c write (iout,*) "distal = ",distal
-c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c & "rescore(",k,j,") =",rescore(k,j)
-c
-c Calculation of sigma from res sim
-c
-c if (odl(k,ii).le.6.0d0) then
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c Other functional forms possible depending on odl(k,ii), eg.
-c
- if (odl(k,ii).le.dist_cut) then
- sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
- else
-#ifdef OLDSIGMA
- sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
- & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
-#else
- sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
- & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
-#endif
-c Following expr replaced by a positive exp argument
-c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
- endif
-c
- sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
-c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
-c
-c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
-c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
- enddo
-c read (ientin,*) sigma_odl(k,ii) ! 1st variant
- enddo
-c lim_odl=ii
-c if (constr_homology.gt.0) call homology_partition
- endif
-c
-c Theta, dihedral and SC retraints
-c
- if (waga_angle.gt.0.0d0) then
-c open (ientin,file=tpl_k_sigma_dih,status='old')
-c do irec=1,maxres-3 ! loop for reading sigma_dih
-c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
-c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
-c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
-c & sigma_dih(k,i+nnt-1)
-c enddo
-c1402 continue
-c close (ientin)
- do i = nnt+3,nct ! right? without parallel.
-c do i=1,nres ! alternative for bounds acc to readpdb?
-c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
- dih(k,i)=phiref(i) ! right?
-c read (ientin,*) sigma_dih(k,i) ! original variant
-c write (iout,*) "dih(",k,i,") =",dih(k,i)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c & "rescore(",k,i-1,") =",rescore(k,i-1),
-c & "rescore(",k,i-2,") =",rescore(k,i-2),
-c & "rescore(",k,i-3,") =",rescore(k,i-3)
-
- sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
- & rescore(k,i-2)+rescore(k,i-3) ! right expression ?
-c
-c write (iout,*) "Raw sigmas for dihedral angle restraints"
-c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
-c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
-c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
-c Instead of res sim other local measure of b/b str reliability possible
- sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
-c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
- if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
-c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
- enddo
- endif
-
- if (waga_theta.gt.0.0d0) then
-c open (ientin,file=tpl_k_sigma_theta,status='old')
-c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
-c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
-c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
-c & sigma_theta(k,i+nnt-1)
-c enddo
-c1403 continue
-c close (ientin)
-
- do i = nnt+2,nct ! right? without parallel.
-c do i = i=1,nres ! alternative for bounds acc to readpdb?
-c do i=ithet_start,ithet_end ! with FG parallel.
- thetatpl(k,i)=thetaref(i)
-c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c & "rescore(",k,i-1,") =",rescore(k,i-1),
-c & "rescore(",k,i-2,") =",rescore(k,i-2)
-c read (ientin,*) sigma_theta(k,i) ! 1st variant
- sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
- & rescore(k,i-2) ! right expression ?
- sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
-
-c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
-c rescore(k,i-2) ! right expression ?
-c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
- if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
- enddo
- endif
-
- if (waga_d.gt.0.0d0) then
-c open (ientin,file=tpl_k_sigma_d,status='old')
-c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
-c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
-c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
-c & sigma_d(k,i+nnt-1)
-c enddo
-c1404 continue
- close (ientin)
-
- do i = nnt,nct ! right? without parallel.
-c do i=2,nres-1 ! alternative for bounds acc to readpdb?
-c do i=loc_start,loc_end ! with FG parallel.
- if (itype(i).eq.10) goto 1 ! right?
- xxtpl(k,i)=xxref(i)
- yytpl(k,i)=yyref(i)
- zztpl(k,i)=zzref(i)
-c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
-c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
-c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
- sigma_d(k,i)=rescore(k,i) ! right expression ?
- sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
-
-c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
-c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
-c read (ientin,*) sigma_d(k,i) ! 1st variant
- if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
- 1 continue
- enddo
- endif
- close(ientin)
- enddo
- if (waga_dist.ne.0.0d0) lim_odl=ii
- if (constr_homology.gt.0) call homology_partition
- if (constr_homology.gt.0) call init_int_table
-cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
-cd & "lim_xx=",lim_xx
-c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
-c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
-c
-c Print restraints
-c
- if (.not.lprn) return
-cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
- if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
- write (iout,*) "Distance restraints from templates"
- do ii=1,lim_odl
- write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
- & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
- enddo
- write (iout,*) "Dihedral angle restraints from templates"
- do i=nnt+3,lim_dih
- write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
- & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
- enddo
- write (iout,*) "Virtual-bond angle restraints from templates"
- do i=nnt+2,lim_theta
- write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
- & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
- enddo
- write (iout,*) "SC restraints from templates"
- do i=nnt,lim_xx
- write(iout,'(i5,10(4f8.2,4x))') i,
- & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
- & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
- enddo
- endif
-c -----------------------------------------------------------------
- return
- end
-
-
->>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit
projects / unres.git / commitdiff