X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=f18e6ebd1ee047393a0f54d52c0f9be1d9fdcc19;hb=d24fd92f8762d3b48a94c699ef90e0ab37d3c204;hp=ab7b8086367cccaf8adacc8113299cc7d1473a27;hpb=ec38ad73a7c0f45390b129510293f6527a5fd8c7;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index ab7b808..f18e6eb 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -98,6 +98,11 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) + write (iout,*) "constr_dist",constr_dist + call readi(controlcard,'NSAXS',nsaxs,0) + call readi(controlcard,'SAXS_MODE',saxs_mode,0) + write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", + & SAXS_MODE call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) call readi(controlcard,'SYM',symetr,1) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours @@ -619,6 +624,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) call reada(weightcard,'WLT',wliptran,0.0D0) + call reada(weightcard,'WSAXS',wsaxs,1.0D0) 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) @@ -642,6 +648,7 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(17)=wbond weights(18)=scal14 weights(21)=wsccor + weights(25)=wsaxs if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, @@ -803,12 +810,16 @@ c print *,'Finished reading pdb data' call gen_side(iti,theta(i+1),alph(i),omeg(i),fail) nsi=nsi+1 enddo +c AL 7/10/16 +c Calculalte only the coordinates of the current sidechain; no need to rebuild +c whole chain + call locate_side_chain(i) if(fail) write(iout,*)'Adding sidechain failed for res ', & i,' after ',nsi,' trials' endif enddo C 10/03/12 Adam: Recalculate coordinates with new side chain positions - call chainbuild +c call chainbuild endif endif if (indpdb.eq.0) then @@ -869,6 +880,7 @@ C 8/13/98 Set limits to generating the dihedral angles phibound(2,i)=pi enddo read (inp,*) ndih_constr + write (iout,*) "ndish_constr",ndih_constr if (ndih_constr.gt.0) then read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) @@ -963,7 +975,7 @@ c---------------------- call MPI_Finalize(MPI_COMM_WORLD,IERROR) stop 'Error reading reference structure' #endif - 39 call chainbuild + 39 call chainbuild_extconf call setup_var czscore call geom_to_var(nvar,coord_exp_zs(1,1)) nstart_sup=nnt @@ -998,6 +1010,8 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb enddo endif endif + write (iout,*) "calling read_saxs_consrtr",nsaxs + if (nsaxs.gt.0) call read_saxs_constr if (constr_homology.gt.0) then @@ -1082,7 +1096,10 @@ C initial geometry. enddo return else + write (iout,*) "Calling read_ang" call read_angles(inp,*36) + write (iout,*) "Calling chainbuild" + call chainbuild_extconf endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' @@ -2117,6 +2134,8 @@ C Get parameter filenames and open the parameter files. open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old', & readonly) @@ -2148,7 +2167,7 @@ c print *,ithep_pdb," opened" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) call getenv_loc('LIPTRANPAR',liptranname) - open (iliptranpar,file=liptranname,status='old',action='read') + open (iliptranpar,file=liptranname,status='old',readonly) #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') @@ -2444,6 +2463,67 @@ CCCC NOW PROPERTIES FOR AFM end c------------------------------------------------------------------------------- + subroutine read_saxs_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + double precision cm(3) +c read(inp,*) nsaxs + write (iout,*) "Calling read_saxs nsaxs",nsaxs + call flush(iout) + if (saxs_mode.eq.0) then +c SAXS distance distribution + do i=1,nsaxs + read(inp,*) distsaxs(i),Psaxs(i) + enddo + Cnorm = 0.0d0 + do i=1,nsaxs + Cnorm = Cnorm + Psaxs(i) + enddo + write (iout,*) "Cnorm",Cnorm + do i=1,nsaxs + Psaxs(i)=Psaxs(i)/Cnorm + enddo + write (iout,*) "Normalized distance distribution from SAXS" + do i=1,nsaxs + write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i) + enddo + else +c SAXS "spheres". + do i=1,nsaxs + read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3) + enddo + do j=1,3 + cm(j)=0.0d0 + enddo + do i=1,nsaxs + do j=1,3 + cm(j)=cm(j)+Csaxs(j,i) + enddo + enddo + do j=1,3 + cm(j)=cm(j)/nsaxs + enddo + do i=1,nsaxs + do j=1,3 + Csaxs(j,i)=Csaxs(j,i)-cm(j) + enddo + enddo + write (iout,*) "SAXS sphere coordinates" + do i=1,nsaxs + write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3) + enddo + endif + return + end +c------------------------------------------------------------------------------- subroutine read_dist_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS'