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 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
99 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
100 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
101 if (lipthick.gt.0.0d0) then
102 bordliptop=(boxzsize+lipthick)/2.0
103 bordlipbot=bordliptop-lipthick
105 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
106 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
107 buflipbot=bordlipbot+lipbufthick
108 bufliptop=bordliptop-lipbufthick
109 if ((lipbufthick*2.0d0).gt.lipthick)
110 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
112 write(iout,*) "bordliptop=",bordliptop
113 write(iout,*) "bordlipbot=",bordlipbot
114 write(iout,*) "bufliptop=",bufliptop
115 write(iout,*) "buflipbot=",buflipbot
116 call readi(controlcard,'SYM',symetr,1)
117 write (iout,*) "DISTCHAINMAX",distchainmax
118 write (iout,*) "delta",delta
119 write (iout,*) "einicheck",einicheck
120 write (iout,*) "rescale_mode",rescale_mode
122 bxfile=index(controlcard,"BXFILE").gt.0
123 cxfile=index(controlcard,"CXFILE").gt.0
124 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
126 histfile=index(controlcard,"HISTFILE").gt.0
127 histout=index(controlcard,"HISTOUT").gt.0
128 entfile=index(controlcard,"ENTFILE").gt.0
129 zscfile=index(controlcard,"ZSCFILE").gt.0
130 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
131 write (iout,*) "with_dihed_constr ",with_dihed_constr
132 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
133 write (iout,*) "with_theta_constr ",with_theta_constr
134 call readi(controlcard,'SHIELD',shield_mode,0)
135 write(iout,*) "shield_mode",shield_mode
137 call readi(controlcard,'TORMODE',tor_mode,0)
138 write(iout,*) "torsional and valence angle mode",tor_mode
139 if (shield_mode.gt.0) then
141 C VSolvSphere the volume of solving sphere
143 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
144 C there will be no distinction between proline peptide group and normal peptide
145 C group in case of shielding parameters
146 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
147 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
148 write (iout,*) VSolvSphere,VSolvSphere_div
149 C long axis of side chain
151 C long_r_sidechain(i)=vbldsc0(1,i)
152 C short_r_sidechain(i)=sigma0(i)
157 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
158 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
159 c if (constr_homology) tole=dmax1(tole,1.5d0)
160 write (iout,*) "with_homology_constr ",with_dihed_constr,
161 & " CONSTR_HOMOLOGY",constr_homology
162 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
163 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
164 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
165 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
166 write (iout,*) "out_template_restr",OUT_TEMPLATE_RESTR
167 dyn_ss=index(controlcard,"DYN_SS").gt.0
168 adaptive = index(controlcard,"ADAPTIVE").gt.0
169 call readi(controlcard,'NSAXS',nsaxs,0)
170 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
171 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
172 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
173 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
174 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
177 c------------------------------------------------------------------------------
178 subroutine read_efree(*)
180 C Read molecular data
184 include 'DIMENSIONS.ZSCOPT'
185 include 'DIMENSIONS.COMPAR'
186 include 'DIMENSIONS.FREE'
187 include 'COMMON.IOUNITS'
188 include 'COMMON.TIME1'
189 include 'COMMON.SBRIDGE'
190 include 'COMMON.CONTROL'
191 include 'COMMON.CHAIN'
192 include 'COMMON.HEADER'
194 include 'COMMON.FREE'
195 character*320 controlcard,ucase
196 integer iparm,ib,i,j,npars
208 call card_concat(controlcard,.true.)
209 call readi(controlcard,'NT',nT_h(iparm),1)
210 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
212 if (nT_h(iparm).gt.MaxT_h) then
213 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
217 replica(iparm)=index(controlcard,"REPLICA").gt.0
218 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
219 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
220 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
221 & replica(iparm)," umbrella ",umbrella(iparm),
222 & " read_iset",read_iset(iparm)
225 call card_concat(controlcard,.true.)
226 call readi(controlcard,'NR',nR(ib,iparm),1)
227 if (umbrella(iparm)) then
230 nRR(ib,iparm)=nR(ib,iparm)
232 if (nR(ib,iparm).gt.MaxR) then
233 write (iout,*) "Error: parameter out of range: NR",
237 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
238 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
239 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
242 call card_concat(controlcard,.true.)
243 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
245 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
250 write (iout,*) "ib",ib," beta_h",
251 & 1.0d0/(0.001987*beta_h(ib,iparm))
252 write (iout,*) "nR",nR(ib,iparm)
253 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
255 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
256 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
268 nR(ib,iparm)=nR(ib,1)
269 if (umbrella(iparm)) then
272 nRR(ib,iparm)=nR(ib,1)
274 beta_h(ib,iparm)=beta_h(ib,1)
276 f(i,ib,iparm)=f(i,ib,1)
278 KH(j,i,ib,iparm)=KH(j,i,ib,1)
279 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
282 replica(iparm)=replica(1)
283 umbrella(iparm)=umbrella(1)
284 read_iset(iparm)=read_iset(1)
292 c-----------------------------------------------------------------------------
293 subroutine read_protein_data(*)
296 include "DIMENSIONS.ZSCOPT"
297 include "DIMENSIONS.FREE"
300 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
303 include "COMMON.CHAIN"
304 include "COMMON.IOUNITS"
305 include "COMMON.PROT"
306 include "COMMON.PROTFILES"
307 include "COMMON.NAMES"
308 include "COMMON.FREE"
309 include "COMMON.OBCINKA"
311 character*16000 controlcard
312 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
322 C Read names of files with conformation data.
323 if (replica(iparm)) then
329 do ii=1,nRR(ib,iparm)
330 write (iout,*) "Parameter set",iparm," temperature",ib,
333 call card_concat(controlcard,.true.)
334 write (iout,*) controlcard(:ilen(controlcard))
335 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
336 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
337 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
338 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
339 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
340 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
341 call reada(controlcard,"TIME_START",
342 & time_start_collect(ii,ib,iparm),0.0d0)
343 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
345 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
346 & " rec_end",rec_end(ii,ib,iparm)
347 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
348 & " time_end",time_end_collect(ii,ib,iparm)
350 if (replica(iparm)) then
351 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
352 write (iout,*) "Number of trajectories",totraj(ii,iparm)
355 if (nfile_bin(ii,ib,iparm).lt.2
356 & .and. nfile_asc(ii,ib,iparm).eq.0
357 & .and. nfile_cx(ii,ib,iparm).eq.0) then
358 write (iout,*) "Error - no action specified!"
361 if (nfile_bin(ii,ib,iparm).gt.0) then
362 call card_concat(controlcard,.false.)
363 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
364 & maxfile_prot,nfile_bin(ii,ib,iparm))
366 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
367 write(iout,*) (protfiles(i,1,ii,ib,iparm),
368 & i=1,nfile_bin(ii,ib,iparm))
371 if (nfile_asc(ii,ib,iparm).gt.0) then
372 call card_concat(controlcard,.false.)
373 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
374 & maxfile_prot,nfile_asc(ii,ib,iparm))
376 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
377 write(iout,*) (protfiles(i,2,ii,ib,iparm),
378 & i=1,nfile_asc(ii,ib,iparm))
380 else if (nfile_cx(ii,ib,iparm).gt.0) then
381 call card_concat(controlcard,.false.)
382 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
383 & maxfile_prot,nfile_cx(ii,ib,iparm))
385 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
386 write(iout,*) (protfiles(i,2,ii,ib,iparm),
387 & i=1,nfile_cx(ii,ib,iparm))
398 c-------------------------------------------------------------------------------
399 subroutine opentmp(islice,iunit,bprotfile_temp)
402 include "DIMENSIONS.ZSCOPT"
403 include "DIMENSIONS.FREE"
406 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
409 include "COMMON.IOUNITS"
410 include "COMMON.PROTFILES"
411 include "COMMON.PROT"
412 include "COMMON.FREE"
413 character*64 bprotfile_temp
414 character*3 liczba,liczba2
421 write (liczba1,'(bz,i2.2)') islice
422 write (liczba,'(bz,i3.3)') me
424 c write (iout,*) "separate_parset ",separate_parset,
426 if (separate_parset) then
427 write (liczba2,'(bz,i3.3)') myparm
428 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
429 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
430 open (iunit,file=bprotfile_temp,status="unknown",
431 & form="unformatted",access="direct",recl=lenrec)
433 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
434 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
435 open (iunit,file=bprotfile_temp,status="unknown",
436 & form="unformatted",access="direct",recl=lenrec)
439 bprotfile_temp = scratchdir(:ilen(scratchdir))//
440 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
441 open (iunit,file=bprotfile_temp,status="unknown",
442 & form="unformatted",access="direct",recl=lenrec)
444 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
449 c-------------------------------------------------------------------------------
450 subroutine read_database(*)
453 include "DIMENSIONS.ZSCOPT"
454 include "DIMENSIONS.FREE"
457 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
460 include "COMMON.CHAIN"
461 include "COMMON.IOUNITS"
462 include "COMMON.PROTFILES"
463 include "COMMON.NAMES"
466 include "COMMON.ENEPS"
467 include "COMMON.PROT"
468 include "COMMON.INTERACT"
469 include "COMMON.FREE"
470 include "COMMON.SBRIDGE"
471 include "COMMON.OBCINKA"
472 real*4 csingle(3,maxres2)
473 character*64 nazwa,bprotfile_temp
476 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
477 & ll(maxslice),mm(maxslice),if
478 integer nrec,nlines,iscor,iunit,islice
479 double precision energ
482 double precision rmsdev,energia(0:max_ene),efree,eini,temp
483 double precision prop(maxQ)
484 integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
485 integer iparm,ib,iib,ir,nprop,nthr,npars
486 double precision etot,time
490 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
491 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
493 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
503 write (iout,*) "nparmset",nparmset
511 if (replica(iparm)) then
518 do iR=1,nRR(ib,iparm)
520 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
526 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
527 c Read conformations from binary DA files (one per batch) and write them to
528 c a binary DA scratchfile.
529 write (liczba,'(bz,i3.3)') me
530 do if=1,nfile_bin(iR,ib,iparm)
531 nazwa=protfiles(if,1,iR,ib,iparm)
532 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
533 open (ientin,file=nazwa,status="old",form="unformatted",
534 & access="direct",recl=lenrec2,err=1111)
537 call opentmp(islice,ientout,bprotfile_temp)
538 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
539 & mm(islice),iR,ib,iparm)
546 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
547 c Read conformations from multiple ASCII int files and write them to a binary
549 do if=1,nfile_asc(iR,ib,iparm)
550 nazwa=protfiles(if,2,iR,ib,iparm)
551 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
552 open(unit=ientin,file=nazwa,status='old',err=1111)
553 write(iout,*) "reading ",nazwa(:ilen(nazwa))
555 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
558 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
559 c Read conformations from cx files and write them to a binary
561 do if=1,nfile_cx(iR,ib,iparm)
562 nazwa=protfiles(if,2,iR,ib,iparm)
563 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
564 write(iout,*) "reading ",nazwa(:ilen(nazwa))
566 print *,"Calling cxread"
567 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
570 write (iout,*) "exit cxread"
576 stot(islice)=stot(islice)+jj(islice)
581 write (iout,*) "IPARM",iparm
584 if (nslice.eq.1) then
586 write (liczba,'(bz,i3.3)') me
587 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
588 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
590 bprotfile_temp = scratchdir(:ilen(scratchdir))//
591 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
593 write(iout,*) mm(1)," conformations read",ll(1),
594 & " conformations written to ",
595 & bprotfile_temp(:ilen(bprotfile_temp))
598 write (liczba1,'(bz,i2.2)') islice
600 write (liczba,'(bz,i3.3)') me
601 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
602 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
604 bprotfile_temp = scratchdir(:ilen(scratchdir))//
605 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
607 write(iout,*) mm(islice)," conformations read",ll(islice),
608 & " conformations written to ",
609 & bprotfile_temp(:ilen(bprotfile_temp))
614 c Check if everyone has the same number of conformations
616 c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
617 c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
619 maxslice_buff=maxslice
621 call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
622 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
627 if (stot(islice).ne.ntot_all(islice,i)) then
628 write (iout,*) "Number of conformations at processor",i,
629 & " differs from that at processor",me,
630 & stot(islice),ntot_all(islice,i)," slice",islice
638 write (iout,*) "Numbers of conformations read by processors"
641 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
643 write (iout,*) "Calculation terminated."
648 ntot(islice)=stot(islice)
652 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
656 c------------------------------------------------------------------------------
657 subroutine card_concat(card,to_upper)
659 include 'DIMENSIONS.ZSCOPT'
660 include "COMMON.IOUNITS"
662 character*80 karta,ucase
666 read (inp,'(a)') karta
667 if (to_upper) karta=ucase(karta)
669 do while (karta(80:80).eq.'&')
670 card=card(:ilen(card)+1)//karta(:79)
671 read (inp,'(a)') karta
672 if (to_upper) karta=ucase(karta)
674 card=card(:ilen(card)+1)//karta
677 c------------------------------------------------------------------------------
678 subroutine readi(rekord,lancuch,wartosc,default)
680 character*(*) rekord,lancuch
681 integer wartosc,default
684 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
689 iread=iread+ilen(lancuch)+1
690 read (rekord(iread:),*) wartosc
693 c----------------------------------------------------------------------------
694 subroutine reada(rekord,lancuch,wartosc,default)
696 character*(*) rekord,lancuch
698 double precision wartosc,default
701 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
706 iread=iread+ilen(lancuch)+1
707 read (rekord(iread:),*) wartosc
710 c----------------------------------------------------------------------------
711 subroutine multreadi(rekord,lancuch,tablica,dim,default)
714 integer tablica(dim),default
715 character*(*) rekord,lancuch
722 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
723 if (iread.eq.0) return
724 iread=iread+ilen(lancuch)+1
725 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
728 c----------------------------------------------------------------------------
729 subroutine multreada(rekord,lancuch,tablica,dim,default)
732 double precision tablica(dim),default
733 character*(*) rekord,lancuch
740 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
741 if (iread.eq.0) return
742 iread=iread+ilen(lancuch)+1
743 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
746 c----------------------------------------------------------------------------
747 subroutine reads(rekord,lancuch,wartosc,default)
749 character*(*) rekord,lancuch,wartosc,default
751 integer ilen,lenlan,lenrec,iread,ireade
757 iread=index(rekord,lancuch(:lenlan)//"=")
758 c print *,"rekord",rekord," lancuch",lancuch
759 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
765 c print *,"iread",iread
766 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
767 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
769 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
771 c print *,"iread",iread
772 if (iread.gt.lenrec) then
777 c print *,"ireade",ireade
778 do while (ireade.lt.lenrec .and.
779 & .not.iblnk(rekord(ireade:ireade)))
782 wartosc=rekord(iread:ireade)
785 c----------------------------------------------------------------------------
786 subroutine multreads(rekord,lancuch,tablica,dim,default)
789 character*(*) rekord,lancuch,tablica(dim),default
791 integer ilen,lenlan,lenrec,iread,ireade
800 iread=index(rekord,lancuch(:lenlan)//"=")
801 c print *,"rekord",rekord," lancuch",lancuch
802 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
803 if (iread.eq.0) return
806 c print *,"iread",iread
807 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
808 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
810 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
812 c print *,"iread",iread
813 if (iread.gt.lenrec) return
815 c print *,"ireade",ireade
816 do while (ireade.lt.lenrec .and.
817 & .not.iblnk(rekord(ireade:ireade)))
820 tablica(i)=rekord(iread:ireade)
824 c----------------------------------------------------------------------------
825 subroutine split_string(rekord,tablica,dim,nsub)
827 integer dim,nsub,i,ii,ll,kk
828 character*(*) tablica(dim)
839 C Find the start of term name
841 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
844 C Parse the name into TABLICA(i) until blank found
845 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
847 tablica(i)(kk:kk)=rekord(ii:ii)
850 if (kk.gt.0) nsub=nsub+1
855 c--------------------------------------------------------------------------------
856 integer function iroof(n,m)
858 if (ii*m .lt. n) ii=ii+1
862 c-------------------------------------------------------------------------------
863 subroutine read_dist_constr
864 implicit real*8 (a-h,o-z)
866 include 'COMMON.CONTROL'
867 include 'COMMON.CHAIN'
868 include 'COMMON.IOUNITS'
869 include 'COMMON.SBRIDGE'
870 include 'COMMON.INTERACT'
871 integer ifrag_(2,100),ipair_(2,100)
872 double precision wfrag_(100),wpair_(100)
873 character*500 controlcard
874 logical normalize,next
876 double precision xlink(4,0:4) /
878 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
879 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
880 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
881 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
882 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
883 write (iout,*) "Calling read_dist_constr"
884 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
886 restr_on_coord=.false.
891 call card_concat(controlcard,.true.)
892 next = index(controlcard,"NEXT").gt.0
893 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
894 write (iout,*) "restr_type",restr_type
895 call readi(controlcard,"NFRAG",nfrag_,0)
896 call readi(controlcard,"NFRAG",nfrag_,0)
897 call readi(controlcard,"NPAIR",npair_,0)
898 call readi(controlcard,"NDIST",ndist_,0)
899 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
900 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
901 if (restr_type.eq.10)
902 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
903 if (restr_type.eq.12)
904 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
905 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
906 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
907 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
908 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
909 normalize = index(controlcard,"NORMALIZE").gt.0
910 write (iout,*) "WBOLTZD",wboltzd
911 write (iout,*) "SCAL_PEAK",scal_peak
912 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
913 write (iout,*) "IFRAG"
915 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
917 write (iout,*) "IPAIR"
919 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
921 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then
923 read(inp,'(a)') pdbfile
925 & "Distance restraints will be constructed from structure ",pdbfile
926 open(ipdbin,file=pdbfile,status='old',err=11)
932 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
933 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
934 & ifrag_(2,i)=nstart_sup+nsup-1
935 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
937 if (wfrag_(i).eq.0.0d0) cycle
938 do j=ifrag_(1,i),ifrag_(2,i)-1
940 c write (iout,*) "j",j," k",k
942 if (restr_type.eq.1) then
948 forcon(nhpb)=wfrag_(i)
949 else if (constr_dist.eq.2) then
950 if (ddjk.le.dist_cut) then
956 forcon(nhpb)=wfrag_(i)
958 else if (restr_type.eq.3) then
964 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
966 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
967 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
972 if (wpair_(i).eq.0.0d0) cycle
980 do j=ifrag_(1,ii),ifrag_(2,ii)
981 do k=ifrag_(1,jj),ifrag_(2,jj)
983 if (restr_type.eq.1) then
989 forcon(nhpb)=wpair_(i)
990 else if (constr_dist.eq.2) then
991 if (ddjk.le.dist_cut) then
997 forcon(nhpb)=wpair_(i)
999 else if (restr_type.eq.3) then
1005 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1007 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
1008 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1014 write (iout,*) "Distance restraints as read from input"
1016 if (restr_type.eq.12) then
1017 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1018 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1019 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1020 & fordepth_peak(nhpb_peak+1),npeak
1021 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
1022 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
1023 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
1024 c & fordepth_peak(nhpb_peak+1),npeak
1025 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
1026 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
1027 nhpb_peak=nhpb_peak+1
1028 irestr_type_peak(nhpb_peak)=12
1029 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
1031 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1032 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
1033 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
1034 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
1035 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
1036 if (ibecarb_peak(nhpb_peak).eq.3) then
1037 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1038 else if (ibecarb_peak(nhpb_peak).eq.2) then
1039 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1040 else if (ibecarb_peak(nhpb_peak).eq.1) then
1041 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
1042 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
1044 else if (restr_type.eq.11) then
1045 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1046 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
1047 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
1048 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
1050 irestr_type(nhpb)=11
1051 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
1052 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1053 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
1054 c if (ibecarb(nhpb).gt.0) then
1055 c ihpb(nhpb)=ihpb(nhpb)+nres
1056 c jhpb(nhpb)=jhpb(nhpb)+nres
1058 if (ibecarb(nhpb).eq.3) then
1059 ihpb(nhpb)=ihpb(nhpb)+nres
1060 else if (ibecarb(nhpb).eq.2) then
1061 ihpb(nhpb)=ihpb(nhpb)+nres
1062 else if (ibecarb(nhpb).eq.1) then
1063 ihpb(nhpb)=ihpb(nhpb)+nres
1064 jhpb(nhpb)=jhpb(nhpb)+nres
1066 else if (restr_type.eq.10) then
1067 c Cross-lonk Markov-like potential
1068 call card_concat(controlcard,.true.)
1069 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
1070 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
1072 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
1073 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
1074 if (index(controlcard,"ZL").gt.0) then
1076 else if (index(controlcard,"ADH").gt.0) then
1078 else if (index(controlcard,"PDH").gt.0) then
1080 else if (index(controlcard,"DSS").gt.0) then
1085 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
1086 & xlink(1,link_type))
1087 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
1088 & xlink(2,link_type))
1089 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
1090 & xlink(3,link_type))
1091 call reada(controlcard,"SIGMA",forcon(nhpb+1),
1092 & xlink(4,link_type))
1093 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
1094 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
1095 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
1096 if (forcon(nhpb+1).le.0.0d0 .or.
1097 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
1099 irestr_type(nhpb)=10
1100 c if (ibecarb(nhpb).gt.0) then
1101 c ihpb(nhpb)=ihpb(nhpb)+nres
1102 c jhpb(nhpb)=jhpb(nhpb)+nres
1104 if (ibecarb(nhpb).eq.3) then
1105 jhpb(nhpb)=jhpb(nhpb)+nres
1106 else if (ibecarb(nhpb).eq.2) then
1107 ihpb(nhpb)=ihpb(nhpb)+nres
1108 else if (ibecarb(nhpb).eq.1) then
1109 ihpb(nhpb)=ihpb(nhpb)+nres
1110 jhpb(nhpb)=jhpb(nhpb)+nres
1112 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1113 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1114 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1118 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
1119 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
1120 if (forcon(nhpb+1).gt.0.0d0) then
1122 if (dhpb1(nhpb).eq.0.0d0) then
1127 c if (ibecarb(nhpb).gt.0) then
1128 c ihpb(nhpb)=ihpb(nhpb)+nres
1129 c jhpb(nhpb)=jhpb(nhpb)+nres
1131 if (ibecarb(nhpb).eq.3) then
1132 ihpb(nhpb)=ihpb(nhpb)+nres
1133 else if (ibecarb(nhpb).eq.2) then
1134 ihpb(nhpb)=ihpb(nhpb)+nres
1135 else if (ibecarb(nhpb).eq.1) then
1136 ihpb(nhpb)=ihpb(nhpb)+nres
1137 jhpb(nhpb)=jhpb(nhpb)+nres
1139 if (dhpb(nhpb).eq.0.0d0)
1140 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1142 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
1143 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
1145 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1146 C if (forcon(nhpb+1).gt.0.0d0) then
1148 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1151 if (restr_type.eq.4) then
1152 write (iout,*) "The BFAC array"
1154 write (iout,'(i5,f10.5)') i,bfac(i)
1157 if (itype(i).eq.ntyp1) cycle
1159 if (itype(j).eq.ntyp1) cycle
1160 if (itype(i).eq.10) then
1165 if (itype(j).eq.10) then
1175 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
1178 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
1179 ihpb(nhpb)=i+nres*ii
1180 jhpb(nhpb)=j+nres*jj
1181 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
1182 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
1183 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
1184 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
1193 if (restr_type.eq.5) then
1194 restr_on_coord=.true.
1196 if (itype(i).eq.ntyp1) cycle
1197 bfac(i)=(scal_bfac/bfac(i))**2
1206 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
1207 & fordepthmax=fordepth(i)
1210 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
1213 if (nhpb.gt.nss) then
1214 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
1215 & "The following",nhpb-nss,
1216 & " distance restraints have been imposed:",
1217 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
1220 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
1221 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
1229 11 write (iout,*)"read_dist_restr: error reading reference structure"