1 subroutine read_general_data(*)
4 include "DIMENSIONS.ZSCOPT"
5 include "DIMENSIONS.FREE"
6 include "COMMON.TORSION"
7 include "COMMON.INTERACT"
8 include "COMMON.IOUNITS"
11 include "COMMON.PROTFILES"
12 include "COMMON.CHAIN"
13 include "COMMON.NAMES"
14 include "COMMON.FFIELD"
15 include "COMMON.ENEPS"
16 include "COMMON.WEIGHTS"
18 include "COMMON.CONTROL"
19 include "COMMON.ENERGIES"
20 include "COMMON.SPLITELE"
21 include "COMMON.SBRIDGE"
22 include "COMMON.SHIELD"
24 character*800 controlcard
25 integer i,j,k,ii,n_ene_found
26 integer ind,itype1,itype2,itypf,itypsc,itypp
33 call card_concat(controlcard,.true.)
34 call readi(controlcard,"N_ENE",n_ene,max_ene)
35 if (n_ene.gt.max_ene) then
36 write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
40 call readi(controlcard,"NPARMSET",nparmset,1)
41 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
42 call readi(controlcard,"IPARMPRINT",iparmprint,1)
43 write (iout,*) "PARMPRINT",iparmprint
44 if (nparmset.gt.max_parm) then
45 write (iout,*) "Error: parameter out of range: NPARMSET",
49 call readi(controlcard,"MAXIT",maxit,5000)
50 call reada(controlcard,"FIMIN",fimin,1.0d-3)
51 call readi(controlcard,"ENSEMBLES",ensembles,0)
52 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
53 write (iout,*) "Number of energy parameter sets",nparmset
54 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
56 if (isampl(i).eq.0) then
57 write (iout,*) "ERROR: isampl is 0 for parmset",i
62 write (iout,*) "MaxSlice",MaxSlice
63 call readi(controlcard,"NSLICE",nslice,1)
65 if (nslice.gt.MaxSlice) then
66 write (iout,*) "Error: parameter out of range: NSLICE",nslice,
70 write (iout,*) "Frequency of storing conformations",
71 & (isampl(i),i=1,nparmset)
72 write (iout,*) "Maxit",maxit," Fimin",fimin
73 call readi(controlcard,"NQ",nQ,1)
75 write (iout,*) "Error: parameter out of range: NQ",nq,
80 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
81 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
82 call reada(controlcard,"DELTA",delta,1.0d-2)
83 call reada(controlcard,"TOLE",tole,1.0d-1)
84 call readi(controlcard,"EINICHECK",einicheck,2)
85 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
86 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
87 call readi(controlcard,"RESCALE",rescale_mode,1)
88 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
89 call readi(controlcard,'TORMODE',tor_mode,0)
90 write(iout,*) "torsional and valence angle mode",tor_mode
91 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
92 call reada(controlcard,'BOXX',boxxsize,100.0d0)
93 call reada(controlcard,'BOXY',boxysize,100.0d0)
94 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
95 c Cutoff range for interactions
96 call reada(controlcard,"R_CUT",r_cut,25.0d0)
97 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
98 write (iout,*) "Cutoff on interactions",r_cut
99 write (iout,*) "lambda",rlamb
100 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
101 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
102 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
103 if (lipthick.gt.0.0d0) then
104 bordliptop=(boxzsize+lipthick)/2.0
105 bordlipbot=bordliptop-lipthick
107 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
108 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
109 buflipbot=bordlipbot+lipbufthick
110 bufliptop=bordliptop-lipbufthick
111 if ((lipbufthick*2.0d0).gt.lipthick)
112 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
114 write(iout,*) "bordliptop=",bordliptop
115 write(iout,*) "bordlipbot=",bordlipbot
116 write(iout,*) "bufliptop=",bufliptop
117 write(iout,*) "buflipbot=",buflipbot
118 call readi(controlcard,'SYM',symetr,1)
119 write (iout,*) "DISTCHAINMAX",distchainmax
120 write (iout,*) "delta",delta
121 write (iout,*) "einicheck",einicheck
122 write (iout,*) "rescale_mode",rescale_mode
124 bxfile=index(controlcard,"BXFILE").gt.0
125 cxfile=index(controlcard,"CXFILE").gt.0
126 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
128 histfile=index(controlcard,"HISTFILE").gt.0
129 histout=index(controlcard,"HISTOUT").gt.0
130 entfile=index(controlcard,"ENTFILE").gt.0
131 zscfile=index(controlcard,"ZSCFILE").gt.0
132 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
133 write (iout,*) "with_dihed_constr ",with_dihed_constr
134 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
135 write (iout,*) "with_theta_constr ",with_theta_constr
136 call readi(controlcard,'SHIELD',shield_mode,0)
137 write(iout,*) "shield_mode",shield_mode
139 call readi(controlcard,'TORMODE',tor_mode,0)
140 write(iout,*) "torsional and valence angle mode",tor_mode
141 if (shield_mode.gt.0) then
143 C VSolvSphere the volume of solving sphere
145 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
146 C there will be no distinction between proline peptide group and normal peptide
147 C group in case of shielding parameters
148 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
149 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
150 write (iout,*) VSolvSphere,VSolvSphere_div
151 C long axis of side chain
153 C long_r_sidechain(i)=vbldsc0(1,i)
154 C short_r_sidechain(i)=sigma0(i)
159 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
160 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
161 c if (constr_homology) tole=dmax1(tole,1.5d0)
162 write (iout,*) "with_homology_constr ",with_dihed_constr,
163 & " CONSTR_HOMOLOGY",constr_homology
164 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
165 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
166 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
167 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
168 write (iout,*) "out_template_restr",OUT_TEMPLATE_RESTR
169 dyn_ss=index(controlcard,"DYN_SS").gt.0
170 adaptive = index(controlcard,"ADAPTIVE").gt.0
171 call readi(controlcard,'NSAXS',nsaxs,0)
172 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
173 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
174 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
175 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
176 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
179 c------------------------------------------------------------------------------
180 subroutine read_efree(*)
182 C Read molecular data
186 include 'DIMENSIONS.ZSCOPT'
187 include 'DIMENSIONS.COMPAR'
188 include 'DIMENSIONS.FREE'
189 include 'COMMON.IOUNITS'
190 include 'COMMON.TIME1'
191 include 'COMMON.SBRIDGE'
192 include 'COMMON.CONTROL'
193 include 'COMMON.CHAIN'
194 include 'COMMON.HEADER'
196 include 'COMMON.FREE'
197 character*320 controlcard,ucase
198 integer iparm,ib,i,j,npars
210 call card_concat(controlcard,.true.)
211 call readi(controlcard,'NT',nT_h(iparm),1)
212 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
214 if (nT_h(iparm).gt.MaxT_h) then
215 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
219 replica(iparm)=index(controlcard,"REPLICA").gt.0
220 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
221 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
222 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
223 & replica(iparm)," umbrella ",umbrella(iparm),
224 & " read_iset",read_iset(iparm)
227 call card_concat(controlcard,.true.)
228 call readi(controlcard,'NR',nR(ib,iparm),1)
229 if (umbrella(iparm)) then
232 nRR(ib,iparm)=nR(ib,iparm)
234 if (nR(ib,iparm).gt.MaxR) then
235 write (iout,*) "Error: parameter out of range: NR",
239 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
240 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
241 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
244 call card_concat(controlcard,.true.)
245 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
247 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
252 write (iout,*) "ib",ib," beta_h",
253 & 1.0d0/(0.001987*beta_h(ib,iparm))
254 write (iout,*) "nR",nR(ib,iparm)
255 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
257 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
258 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
270 nR(ib,iparm)=nR(ib,1)
271 if (umbrella(iparm)) then
274 nRR(ib,iparm)=nR(ib,1)
276 beta_h(ib,iparm)=beta_h(ib,1)
278 f(i,ib,iparm)=f(i,ib,1)
280 KH(j,i,ib,iparm)=KH(j,i,ib,1)
281 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
284 replica(iparm)=replica(1)
285 umbrella(iparm)=umbrella(1)
286 read_iset(iparm)=read_iset(1)
294 c-----------------------------------------------------------------------------
295 subroutine read_protein_data(*)
298 include "DIMENSIONS.ZSCOPT"
299 include "DIMENSIONS.FREE"
302 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
305 include "COMMON.CHAIN"
306 include "COMMON.IOUNITS"
307 include "COMMON.PROT"
308 include "COMMON.PROTFILES"
309 include "COMMON.NAMES"
310 include "COMMON.FREE"
311 include "COMMON.OBCINKA"
313 character*16000 controlcard
314 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
324 C Read names of files with conformation data.
325 if (replica(iparm)) then
331 do ii=1,nRR(ib,iparm)
332 write (iout,*) "Parameter set",iparm," temperature",ib,
335 call card_concat(controlcard,.true.)
336 write (iout,*) controlcard(:ilen(controlcard))
337 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
338 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
339 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
340 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
341 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
342 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
343 call reada(controlcard,"TIME_START",
344 & time_start_collect(ii,ib,iparm),0.0d0)
345 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
347 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
348 & " rec_end",rec_end(ii,ib,iparm)
349 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
350 & " time_end",time_end_collect(ii,ib,iparm)
352 if (replica(iparm)) then
353 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
354 write (iout,*) "Number of trajectories",totraj(ii,iparm)
357 if (nfile_bin(ii,ib,iparm).lt.2
358 & .and. nfile_asc(ii,ib,iparm).eq.0
359 & .and. nfile_cx(ii,ib,iparm).eq.0) then
360 write (iout,*) "Error - no action specified!"
363 if (nfile_bin(ii,ib,iparm).gt.0) then
364 call card_concat(controlcard,.false.)
365 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
366 & maxfile_prot,nfile_bin(ii,ib,iparm))
368 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
369 write(iout,*) (protfiles(i,1,ii,ib,iparm),
370 & i=1,nfile_bin(ii,ib,iparm))
373 if (nfile_asc(ii,ib,iparm).gt.0) then
374 call card_concat(controlcard,.false.)
375 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
376 & maxfile_prot,nfile_asc(ii,ib,iparm))
378 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
379 write(iout,*) (protfiles(i,2,ii,ib,iparm),
380 & i=1,nfile_asc(ii,ib,iparm))
382 else if (nfile_cx(ii,ib,iparm).gt.0) then
383 call card_concat(controlcard,.false.)
384 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
385 & maxfile_prot,nfile_cx(ii,ib,iparm))
387 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
388 write(iout,*) (protfiles(i,2,ii,ib,iparm),
389 & i=1,nfile_cx(ii,ib,iparm))
400 c-------------------------------------------------------------------------------
401 subroutine opentmp(islice,iunit,bprotfile_temp)
404 include "DIMENSIONS.ZSCOPT"
405 include "DIMENSIONS.FREE"
408 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
411 include "COMMON.IOUNITS"
412 include "COMMON.PROTFILES"
413 include "COMMON.PROT"
414 include "COMMON.FREE"
415 character*64 bprotfile_temp
416 character*3 liczba,liczba2
423 write (liczba1,'(bz,i2.2)') islice
424 write (liczba,'(bz,i3.3)') me
426 c write (iout,*) "separate_parset ",separate_parset,
428 if (separate_parset) then
429 write (liczba2,'(bz,i3.3)') myparm
430 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
431 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
432 open (iunit,file=bprotfile_temp,status="unknown",
433 & form="unformatted",access="direct",recl=lenrec)
435 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
436 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
437 open (iunit,file=bprotfile_temp,status="unknown",
438 & form="unformatted",access="direct",recl=lenrec)
441 bprotfile_temp = scratchdir(:ilen(scratchdir))//
442 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
443 open (iunit,file=bprotfile_temp,status="unknown",
444 & form="unformatted",access="direct",recl=lenrec)
446 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
451 c-------------------------------------------------------------------------------
452 subroutine read_database(*)
455 include "DIMENSIONS.ZSCOPT"
456 include "DIMENSIONS.FREE"
459 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
462 include "COMMON.CHAIN"
463 include "COMMON.IOUNITS"
464 include "COMMON.PROTFILES"
465 include "COMMON.NAMES"
468 include "COMMON.ENEPS"
469 include "COMMON.PROT"
470 include "COMMON.INTERACT"
471 include "COMMON.FREE"
472 include "COMMON.SBRIDGE"
473 include "COMMON.OBCINKA"
474 real*4 csingle(3,maxres2)
475 character*64 nazwa,bprotfile_temp
478 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
479 & ll(maxslice),mm(maxslice),if
480 integer nrec,nlines,iscor,iunit,islice
481 double precision energ
484 double precision rmsdev,energia(0:max_ene),efree,eini,temp
485 double precision prop(maxQ)
486 integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
487 integer iparm,ib,iib,ir,nprop,nthr,npars
488 double precision etot,time
492 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
493 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
495 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
505 write (iout,*) "nparmset",nparmset
513 if (replica(iparm)) then
520 do iR=1,nRR(ib,iparm)
522 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
528 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
529 c Read conformations from binary DA files (one per batch) and write them to
530 c a binary DA scratchfile.
531 write (liczba,'(bz,i3.3)') me
532 do if=1,nfile_bin(iR,ib,iparm)
533 nazwa=protfiles(if,1,iR,ib,iparm)
534 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
535 open (ientin,file=nazwa,status="old",form="unformatted",
536 & access="direct",recl=lenrec2,err=1111)
539 call opentmp(islice,ientout,bprotfile_temp)
540 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
541 & mm(islice),iR,ib,iparm)
548 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
549 c Read conformations from multiple ASCII int files and write them to a binary
551 do if=1,nfile_asc(iR,ib,iparm)
552 nazwa=protfiles(if,2,iR,ib,iparm)
553 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
554 open(unit=ientin,file=nazwa,status='old',err=1111)
555 write(iout,*) "reading ",nazwa(:ilen(nazwa))
557 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
560 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
561 c Read conformations from cx files and write them to a binary
563 do if=1,nfile_cx(iR,ib,iparm)
564 nazwa=protfiles(if,2,iR,ib,iparm)
565 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
566 write(iout,*) "reading ",nazwa(:ilen(nazwa))
568 print *,"Calling cxread"
569 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
572 write (iout,*) "exit cxread"
578 stot(islice)=stot(islice)+jj(islice)
583 write (iout,*) "IPARM",iparm
586 if (nslice.eq.1) then
588 write (liczba,'(bz,i3.3)') me
589 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
590 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
592 bprotfile_temp = scratchdir(:ilen(scratchdir))//
593 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
595 write(iout,*) mm(1)," conformations read",ll(1),
596 & " conformations written to ",
597 & bprotfile_temp(:ilen(bprotfile_temp))
600 write (liczba1,'(bz,i2.2)') islice
602 write (liczba,'(bz,i3.3)') me
603 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
604 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
606 bprotfile_temp = scratchdir(:ilen(scratchdir))//
607 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
609 write(iout,*) mm(islice)," conformations read",ll(islice),
610 & " conformations written to ",
611 & bprotfile_temp(:ilen(bprotfile_temp))
616 c Check if everyone has the same number of conformations
618 c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
619 c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
621 maxslice_buff=maxslice
623 call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
624 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
629 if (stot(islice).ne.ntot_all(islice,i)) then
630 write (iout,*) "Number of conformations at processor",i,
631 & " differs from that at processor",me,
632 & stot(islice),ntot_all(islice,i)," slice",islice
640 write (iout,*) "Numbers of conformations read by processors"
643 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
645 write (iout,*) "Calculation terminated."
650 ntot(islice)=stot(islice)
654 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
658 c------------------------------------------------------------------------------
659 subroutine card_concat(card,to_upper)
661 include 'DIMENSIONS.ZSCOPT'
662 include "COMMON.IOUNITS"
664 character*80 karta,ucase
668 read (inp,'(a)') karta
669 if (to_upper) karta=ucase(karta)
671 do while (karta(80:80).eq.'&')
672 card=card(:ilen(card)+1)//karta(:79)
673 read (inp,'(a)') karta
674 if (to_upper) karta=ucase(karta)
676 card=card(:ilen(card)+1)//karta
679 c------------------------------------------------------------------------------
680 subroutine readi(rekord,lancuch,wartosc,default)
682 character*(*) rekord,lancuch
683 integer wartosc,default
686 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
691 iread=iread+ilen(lancuch)+1
692 read (rekord(iread:),*) wartosc
695 c----------------------------------------------------------------------------
696 subroutine reada(rekord,lancuch,wartosc,default)
698 character*(*) rekord,lancuch
700 double precision wartosc,default
703 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
708 iread=iread+ilen(lancuch)+1
709 read (rekord(iread:),*) wartosc
712 c----------------------------------------------------------------------------
713 subroutine multreadi(rekord,lancuch,tablica,dim,default)
716 integer tablica(dim),default
717 character*(*) rekord,lancuch
724 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
725 if (iread.eq.0) return
726 iread=iread+ilen(lancuch)+1
727 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
730 c----------------------------------------------------------------------------
731 subroutine multreada(rekord,lancuch,tablica,dim,default)
734 double precision tablica(dim),default
735 character*(*) rekord,lancuch
742 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
743 if (iread.eq.0) return
744 iread=iread+ilen(lancuch)+1
745 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
748 c----------------------------------------------------------------------------
749 subroutine reads(rekord,lancuch,wartosc,default)
751 character*(*) rekord,lancuch,wartosc,default
753 integer ilen,lenlan,lenrec,iread,ireade
759 iread=index(rekord,lancuch(:lenlan)//"=")
760 c print *,"rekord",rekord," lancuch",lancuch
761 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
767 c print *,"iread",iread
768 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
769 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
771 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
773 c print *,"iread",iread
774 if (iread.gt.lenrec) then
779 c print *,"ireade",ireade
780 do while (ireade.lt.lenrec .and.
781 & .not.iblnk(rekord(ireade:ireade)))
784 wartosc=rekord(iread:ireade)
787 c----------------------------------------------------------------------------
788 subroutine multreads(rekord,lancuch,tablica,dim,default)
791 character*(*) rekord,lancuch,tablica(dim),default
793 integer ilen,lenlan,lenrec,iread,ireade
802 iread=index(rekord,lancuch(:lenlan)//"=")
803 c print *,"rekord",rekord," lancuch",lancuch
804 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
805 if (iread.eq.0) return
808 c print *,"iread",iread
809 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
810 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
812 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
814 c print *,"iread",iread
815 if (iread.gt.lenrec) return
817 c print *,"ireade",ireade
818 do while (ireade.lt.lenrec .and.
819 & .not.iblnk(rekord(ireade:ireade)))
822 tablica(i)=rekord(iread:ireade)
826 c----------------------------------------------------------------------------
827 subroutine split_string(rekord,tablica,dim,nsub)
829 integer dim,nsub,i,ii,ll,kk
830 character*(*) tablica(dim)
841 C Find the start of term name
843 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
846 C Parse the name into TABLICA(i) until blank found
847 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
849 tablica(i)(kk:kk)=rekord(ii:ii)
852 if (kk.gt.0) nsub=nsub+1
857 c--------------------------------------------------------------------------------
858 integer function iroof(n,m)
860 if (ii*m .lt. n) ii=ii+1
864 c-------------------------------------------------------------------------------
865 subroutine read_dist_constr
866 implicit real*8 (a-h,o-z)
868 include 'COMMON.CONTROL'
869 include 'COMMON.CHAIN'
870 include 'COMMON.IOUNITS'
871 include 'COMMON.SBRIDGE'
872 include 'COMMON.INTERACT'
873 integer ifrag_(2,100),ipair_(2,100)
874 double precision wfrag_(100),wpair_(100)
875 character*500 controlcard
876 logical normalize,next
878 double precision xlink(4,0:4) /
880 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
881 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
882 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
883 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
884 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
885 write (iout,*) "Calling read_dist_constr"
886 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
888 restr_on_coord=.false.
893 call card_concat(controlcard,.true.)
894 next = index(controlcard,"NEXT").gt.0
895 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
896 write (iout,*) "restr_type",restr_type
897 call readi(controlcard,"NFRAG",nfrag_,0)
898 call readi(controlcard,"NFRAG",nfrag_,0)
899 call readi(controlcard,"NPAIR",npair_,0)
900 call readi(controlcard,"NDIST",ndist_,0)
901 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
902 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
903 if (restr_type.eq.10)
904 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
905 if (restr_type.eq.12)
906 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
907 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
908 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
909 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
910 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
911 normalize = index(controlcard,"NORMALIZE").gt.0
912 write (iout,*) "WBOLTZD",wboltzd
913 write (iout,*) "SCAL_PEAK",scal_peak
914 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
915 write (iout,*) "IFRAG"
917 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
919 write (iout,*) "IPAIR"
921 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
923 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
925 read(inp,'(a)') pdbfile
927 & "Distance restraints will be constructed from structure ",pdbfile
928 open(ipdbin,file=pdbfile,status='old',err=11)
934 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
935 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
936 & ifrag_(2,i)=nstart_sup+nsup-1
937 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
939 if (wfrag_(i).eq.0.0d0) cycle
940 do j=ifrag_(1,i),ifrag_(2,i)-1
942 c write (iout,*) "j",j," k",k
944 if (restr_type.eq.1) then
950 forcon(nhpb)=wfrag_(i)
951 else if (constr_dist.eq.2) then
952 if (ddjk.le.dist_cut) then
958 forcon(nhpb)=wfrag_(i)
960 else if (restr_type.eq.3) then
966 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
968 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
969 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
974 if (wpair_(i).eq.0.0d0) cycle
982 do j=ifrag_(1,ii),ifrag_(2,ii)
983 do k=ifrag_(1,jj),ifrag_(2,jj)
985 if (restr_type.eq.1) then
991 forcon(nhpb)=wpair_(i)
992 else if (constr_dist.eq.2) then
993 if (ddjk.le.dist_cut) then
999 forcon(nhpb)=wpair_(i)
1001 else if (restr_type.eq.3) then
1007 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1009 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1010 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1016 write (iout,*) "Distance restraints as read from input"
1018 if (restr_type.eq.12) then
1019 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1020 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1021 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1022 & fordepth_peak(nhpb_peak+1),npeak
1023 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1024 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1025 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1026 c & fordepth_peak(nhpb_peak+1),npeak
1027 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1028 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1029 nhpb_peak=nhpb_peak+1
1030 irestr_type_peak(nhpb_peak)=12
1031 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1033 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1034 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1035 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1036 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1037 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1038 if (ibecarb_peak(nhpb_peak).eq.3) then
1039 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1040 else if (ibecarb_peak(nhpb_peak).eq.2) then
1041 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1042 else if (ibecarb_peak(nhpb_peak).eq.1) then
1043 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1044 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1046 else if (restr_type.eq.11) then
1047 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1048 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1049 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1050 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1052 irestr_type(nhpb)=11
1053 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1054 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1055 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1056 c if (ibecarb(nhpb).gt.0) then
1057 c ihpb(nhpb)=ihpb(nhpb)+nres
1058 c jhpb(nhpb)=jhpb(nhpb)+nres
1060 if (ibecarb(nhpb).eq.3) then
1061 ihpb(nhpb)=ihpb(nhpb)+nres
1062 else if (ibecarb(nhpb).eq.2) then
1063 ihpb(nhpb)=ihpb(nhpb)+nres
1064 else if (ibecarb(nhpb).eq.1) then
1065 ihpb(nhpb)=ihpb(nhpb)+nres
1066 jhpb(nhpb)=jhpb(nhpb)+nres
1068 else if (restr_type.eq.10) then
1069 c Cross-lonk Markov-like potential
1070 call card_concat(controlcard,.true.)
1071 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1072 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1074 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1075 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1076 if (index(controlcard,"ZL").gt.0) then
1078 else if (index(controlcard,"ADH").gt.0) then
1080 else if (index(controlcard,"PDH").gt.0) then
1082 else if (index(controlcard,"DSS").gt.0) then
1087 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1088 & xlink(1,link_type))
1089 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1090 & xlink(2,link_type))
1091 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1092 & xlink(3,link_type))
1093 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1094 & xlink(4,link_type))
1095 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1096 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1097 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1098 if (forcon(nhpb+1).le.0.0d0 .or.
1099 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1101 irestr_type(nhpb)=10
1102 c if (ibecarb(nhpb).gt.0) then
1103 c ihpb(nhpb)=ihpb(nhpb)+nres
1104 c jhpb(nhpb)=jhpb(nhpb)+nres
1106 if (ibecarb(nhpb).eq.3) then
1107 jhpb(nhpb)=jhpb(nhpb)+nres
1108 else if (ibecarb(nhpb).eq.2) then
1109 ihpb(nhpb)=ihpb(nhpb)+nres
1110 else if (ibecarb(nhpb).eq.1) then
1111 ihpb(nhpb)=ihpb(nhpb)+nres
1112 jhpb(nhpb)=jhpb(nhpb)+nres
1114 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1115 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1116 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1120 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1121 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1122 if (forcon(nhpb+1).gt.0.0d0) then
1124 if (dhpb1(nhpb).eq.0.0d0) then
1129 c if (ibecarb(nhpb).gt.0) then
1130 c ihpb(nhpb)=ihpb(nhpb)+nres
1131 c jhpb(nhpb)=jhpb(nhpb)+nres
1133 if (ibecarb(nhpb).eq.3) then
1134 ihpb(nhpb)=ihpb(nhpb)+nres
1135 else if (ibecarb(nhpb).eq.2) then
1136 ihpb(nhpb)=ihpb(nhpb)+nres
1137 else if (ibecarb(nhpb).eq.1) then
1138 ihpb(nhpb)=ihpb(nhpb)+nres
1139 jhpb(nhpb)=jhpb(nhpb)+nres
1141 if (dhpb(nhpb).eq.0.0d0)
1142 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1144 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1145 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1147 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1148 C if (forcon(nhpb+1).gt.0.0d0) then
1150 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1153 if (restr_type.eq.4) then
1154 write (iout,*) "The BFAC array"
1156 write (iout,'(i5,f10.5)') i,bfac(i)
1159 if (itype(i).eq.ntyp1) cycle
1161 if (itype(j).eq.ntyp1) cycle
1162 if (itype(i).eq.10) then
1167 if (itype(j).eq.10) then
1177 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1180 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1181 ihpb(nhpb)=i+nres*ii
1182 jhpb(nhpb)=j+nres*jj
1183 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1184 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1185 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1186 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1195 if (restr_type.eq.5) then
1196 restr_on_coord=.true.
1198 if (itype(i).eq.ntyp1) cycle
1199 bfac(i)=(scal_bfac/bfac(i))**2
1208 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1209 & fordepthmax=fordepth(i)
1212 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1215 if (nhpb.gt.nss) then
1216 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1217 & "The following",nhpb-nss,
1218 & " distance restraints have been imposed:",
1219 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1222 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1223 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1231 11 write (iout,*)"read_dist_restr: error reading reference structure"