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 character*800 controlcard
21 integer i,j,k,ii,n_ene_found
22 integer ind,itype1,itype2,itypf,itypsc,itypp
29 call card_concat(controlcard,.true.)
30 call readi(controlcard,"N_ENE",n_ene,max_ene)
31 if (n_ene.gt.max_ene) then
32 write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
36 call readi(controlcard,"NPARMSET",nparmset,1)
37 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
38 call readi(controlcard,"IPARMPRINT",iparmprint,1)
39 write (iout,*) "PARMPRINT",iparmprint
40 if (nparmset.gt.max_parm) then
41 write (iout,*) "Error: parameter out of range: NPARMSET",
47 energy_dec=index(controlcard,"ENERGY_DEC").gt.0
48 >>>>>>> e183793... Added src_MD-M-newcorr (Adasko's source) and src-NEWSC of WHAM (with Momo's SCSC potentials)
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)
55 write (iout,*) "MaxSlice",MaxSlice
56 call readi(controlcard,"NSLICE",nslice,1)
58 if (nslice.gt.MaxSlice) then
59 write (iout,*) "Error: parameter out of range: NSLICE",nslice,
63 write (iout,*) "Frequency of storing conformations",
64 & (isampl(i),i=1,nparmset)
65 write (iout,*) "Maxit",maxit," Fimin",fimin
66 call readi(controlcard,"NQ",nQ,1)
68 write (iout,*) "Error: parameter out of range: NQ",nq,
73 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
74 call reada(controlcard,"DELTA",delta,1.0d-2)
75 call readi(controlcard,"EINICHECK",einicheck,2)
76 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
77 call readi(controlcard,"NGRIDT",NGridT,400)
78 call reada(controlcard,"STARTGRIDT",StartGridT,200.0d0)
79 call reada(controlcard,"DELTA_T",Delta_T,1.0d0)
80 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
81 call readi(controlcard,"RESCALE",rescale_mode,1)
82 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
83 write (iout,*) "delta",delta
84 write (iout,*) "einicheck",einicheck
85 write (iout,*) "rescale_mode",rescale_mode
87 bxfile=index(controlcard,"BXFILE").gt.0
88 cxfile=index(controlcard,"CXFILE").gt.0
89 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
91 histfile=index(controlcard,"HISTFILE").gt.0
92 histout=index(controlcard,"HISTOUT").gt.0
93 entfile=index(controlcard,"ENTFILE").gt.0
94 zscfile=index(controlcard,"ZSCFILE").gt.0
95 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
96 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
97 write (iout,*) "with_dihed_constr ",with_dihed_constr,
98 & " CONSTR_DIST",constr_dist
99 refstr = index(controlcard,'REFSTR').gt.0
100 pdbref = index(controlcard,'PDBREF').gt.0
104 c------------------------------------------------------------------------------
105 subroutine read_efree(*)
107 C Read molecular data
111 include 'DIMENSIONS.ZSCOPT'
112 include 'DIMENSIONS.COMPAR'
113 include 'DIMENSIONS.FREE'
114 include 'COMMON.IOUNITS'
115 include 'COMMON.TIME1'
116 include 'COMMON.SBRIDGE'
117 include 'COMMON.CONTROL'
118 include 'COMMON.CHAIN'
119 include 'COMMON.HEADER'
121 include 'COMMON.FREE'
122 character*320 controlcard,ucase
123 integer iparm,ib,i,j,npars
135 call card_concat(controlcard,.true.)
136 call readi(controlcard,'NT',nT_h(iparm),1)
137 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
139 if (nT_h(iparm).gt.MaxT_h) then
140 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
144 replica(iparm)=index(controlcard,"REPLICA").gt.0
145 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
146 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
147 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
148 & replica(iparm)," umbrella ",umbrella(iparm),
149 & " read_iset",read_iset(iparm)
152 call card_concat(controlcard,.true.)
153 call readi(controlcard,'NR',nR(ib,iparm),1)
154 if (umbrella(iparm)) then
157 nRR(ib,iparm)=nR(ib,iparm)
159 if (nR(ib,iparm).gt.MaxR) then
160 write (iout,*) "Error: parameter out of range: NR",
164 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
165 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
166 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
169 call card_concat(controlcard,.true.)
170 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
172 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
177 write (iout,*) "ib",ib," beta_h",
178 & 1.0d0/(0.001987*beta_h(ib,iparm))
179 write (iout,*) "nR",nR(ib,iparm)
180 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
182 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
183 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
195 nR(ib,iparm)=nR(ib,1)
196 if (umbrella(iparm)) then
199 nRR(ib,iparm)=nR(ib,1)
201 beta_h(ib,iparm)=beta_h(ib,1)
203 f(i,ib,iparm)=f(i,ib,1)
205 KH(j,i,ib,iparm)=KH(j,i,ib,1)
206 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
209 replica(iparm)=replica(1)
210 umbrella(iparm)=umbrella(1)
211 read_iset(iparm)=read_iset(1)
219 c-----------------------------------------------------------------------------
220 subroutine read_protein_data(*)
223 include "DIMENSIONS.ZSCOPT"
224 include "DIMENSIONS.FREE"
227 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
230 include "COMMON.CHAIN"
231 include "COMMON.IOUNITS"
232 include "COMMON.PROT"
233 include "COMMON.PROTFILES"
234 include "COMMON.NAMES"
235 include "COMMON.FREE"
236 include "COMMON.OBCINKA"
238 character*16000 controlcard
239 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
249 C Read names of files with conformation data.
250 if (replica(iparm)) then
256 do ii=1,nRR(ib,iparm)
257 write (iout,*) "Parameter set",iparm," temperature",ib,
260 call card_concat(controlcard,.true.)
261 write (iout,*) controlcard(:ilen(controlcard))
262 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
263 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
264 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
265 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
266 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
267 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
268 call reada(controlcard,"TIME_START",
269 & time_start_collect(ii,ib,iparm),0.0d0)
270 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
272 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
273 & " rec_end",rec_end(ii,ib,iparm)
274 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
275 & " time_end",time_end_collect(ii,ib,iparm)
277 if (replica(iparm)) then
278 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
279 write (iout,*) "Number of trajectories",totraj(ii,iparm)
282 if (nfile_bin(ii,ib,iparm).lt.2
283 & .and. nfile_asc(ii,ib,iparm).eq.0
284 & .and. nfile_cx(ii,ib,iparm).eq.0) then
285 write (iout,*) "Error - no action specified!"
288 if (nfile_bin(ii,ib,iparm).gt.0) then
289 call card_concat(controlcard,.false.)
290 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
291 & maxfile_prot,nfile_bin(ii,ib,iparm))
293 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
294 write(iout,*) (protfiles(i,1,ii,ib,iparm),
295 & i=1,nfile_bin(ii,ib,iparm))
298 if (nfile_asc(ii,ib,iparm).gt.0) then
299 call card_concat(controlcard,.false.)
300 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
301 & maxfile_prot,nfile_asc(ii,ib,iparm))
303 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
304 write(iout,*) (protfiles(i,2,ii,ib,iparm),
305 & i=1,nfile_asc(ii,ib,iparm))
307 else if (nfile_cx(ii,ib,iparm).gt.0) then
308 call card_concat(controlcard,.false.)
309 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
310 & maxfile_prot,nfile_cx(ii,ib,iparm))
312 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
313 write(iout,*) (protfiles(i,2,ii,ib,iparm),
314 & i=1,nfile_cx(ii,ib,iparm))
325 c-------------------------------------------------------------------------------
326 subroutine opentmp(islice,iunit,bprotfile_temp)
329 include "DIMENSIONS.ZSCOPT"
330 include "DIMENSIONS.FREE"
333 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
336 include "COMMON.IOUNITS"
337 include "COMMON.PROTFILES"
338 include "COMMON.PROT"
339 include "COMMON.FREE"
340 character*64 bprotfile_temp
341 character*3 liczba,liczba2
348 write (liczba1,'(bz,i2.2)') islice
349 write (liczba,'(bz,i3.3)') me
351 c write (iout,*) "separate_parset ",separate_parset,
353 if (separate_parset) then
354 write (liczba2,'(bz,i3.3)') myparm
355 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
356 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
357 open (iunit,file=bprotfile_temp,status="unknown",
358 & form="unformatted",access="direct",recl=lenrec)
360 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
361 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
362 open (iunit,file=bprotfile_temp,status="unknown",
363 & form="unformatted",access="direct",recl=lenrec)
366 bprotfile_temp = scratchdir(:ilen(scratchdir))//
367 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
368 open (iunit,file=bprotfile_temp,status="unknown",
369 & form="unformatted",access="direct",recl=lenrec)
371 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
376 c-------------------------------------------------------------------------------
377 subroutine read_database(*)
380 include "DIMENSIONS.ZSCOPT"
381 include "DIMENSIONS.FREE"
384 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
387 include "COMMON.CHAIN"
388 include "COMMON.IOUNITS"
389 include "COMMON.PROTFILES"
390 include "COMMON.NAMES"
393 include "COMMON.ENEPS"
394 include "COMMON.PROT"
395 include "COMMON.INTERACT"
396 include "COMMON.FREE"
397 include "COMMON.SBRIDGE"
398 include "COMMON.OBCINKA"
399 real*4 csingle(3,maxres2)
400 character*64 nazwa,bprotfile_temp
403 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
404 & ll(maxslice),mm(maxslice),if
405 integer nrec,nlines,iscor,iunit,islice
406 double precision energ
409 double precision rmsdev,energia(0:max_ene),efree,eini,temp
410 double precision prop(maxQ)
411 integer ntot_all(maxslice,0:maxprocs-1)
412 integer iparm,ib,iib,ir,nprop,nthr,npars
413 double precision etot,time
417 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
418 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
420 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
430 write (iout,*) "nparmset",nparmset
438 if (replica(iparm)) then
445 do iR=1,nRR(ib,iparm)
447 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
453 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
454 c Read conformations from binary DA files (one per batch) and write them to
455 c a binary DA scratchfile.
456 write (liczba,'(bz,i3.3)') me
457 do if=1,nfile_bin(iR,ib,iparm)
458 nazwa=protfiles(if,1,iR,ib,iparm)
459 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
460 open (ientin,file=nazwa,status="old",form="unformatted",
461 & access="direct",recl=lenrec2,err=1111)
464 call opentmp(islice,ientout,bprotfile_temp)
465 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
466 & mm(islice),iR,ib,iparm)
473 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
474 c Read conformations from multiple ASCII int files and write them to a binary
476 do if=1,nfile_asc(iR,ib,iparm)
477 nazwa=protfiles(if,2,iR,ib,iparm)
478 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
479 open(unit=ientin,file=nazwa,status='old',err=1111)
480 write(iout,*) "reading ",nazwa(:ilen(nazwa))
482 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
485 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
486 c Read conformations from cx files and write them to a binary
488 do if=1,nfile_cx(iR,ib,iparm)
489 nazwa=protfiles(if,2,iR,ib,iparm)
490 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
491 write(iout,*) "reading ",nazwa(:ilen(nazwa))
493 print *,"Calling cxread"
494 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
497 write (iout,*) "exit cxread"
503 stot(islice)=stot(islice)+jj(islice)
508 write (iout,*) "IPARM",iparm
511 if (nslice.eq.1) then
513 write (liczba,'(bz,i3.3)') me
514 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
515 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
517 bprotfile_temp = scratchdir(:ilen(scratchdir))//
518 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
520 write(iout,*) mm(1)," conformations read",ll(1),
521 & " conformations written to ",
522 & bprotfile_temp(:ilen(bprotfile_temp))
525 write (liczba1,'(bz,i2.2)') islice
527 write (liczba,'(bz,i3.3)') me
528 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
529 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
531 bprotfile_temp = scratchdir(:ilen(scratchdir))//
532 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
534 write(iout,*) mm(islice)," conformations read",ll(islice),
535 & " conformations written to ",
536 & bprotfile_temp(:ilen(bprotfile_temp))
541 c Check if everyone has the same number of conformations
542 call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
543 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
548 if (stot(islice).ne.ntot_all(islice,i)) then
549 write (iout,*) "Number of conformations at processor",i,
550 & " differs from that at processor",me,
551 & stot(islice),ntot_all(islice,i)," slice",islice
559 write (iout,*) "Numbers of conformations read by processors"
562 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
564 write (iout,*) "Calculation terminated."
569 ntot(islice)=stot(islice)
573 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
577 c------------------------------------------------------------------------------
578 subroutine card_concat(card,to_upper)
580 include 'DIMENSIONS.ZSCOPT'
581 include "COMMON.IOUNITS"
583 character*80 karta,ucase
587 read (inp,'(a)') karta
588 if (to_upper) karta=ucase(karta)
590 do while (karta(80:80).eq.'&')
591 card=card(:ilen(card)+1)//karta(:79)
592 read (inp,'(a)') karta
593 if (to_upper) karta=ucase(karta)
595 card=card(:ilen(card)+1)//karta
598 c------------------------------------------------------------------------------
599 subroutine readi(rekord,lancuch,wartosc,default)
601 character*(*) rekord,lancuch
602 integer wartosc,default
605 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
610 iread=iread+ilen(lancuch)+1
611 read (rekord(iread:),*) wartosc
614 c----------------------------------------------------------------------------
615 subroutine reada(rekord,lancuch,wartosc,default)
617 character*(*) rekord,lancuch
619 double precision wartosc,default
622 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
627 iread=iread+ilen(lancuch)+1
628 read (rekord(iread:),*) wartosc
631 c----------------------------------------------------------------------------
632 subroutine multreadi(rekord,lancuch,tablica,dim,default)
635 integer tablica(dim),default
636 character*(*) rekord,lancuch
643 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
644 if (iread.eq.0) return
645 iread=iread+ilen(lancuch)+1
646 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
649 c----------------------------------------------------------------------------
650 subroutine multreada(rekord,lancuch,tablica,dim,default)
653 double precision tablica(dim),default
654 character*(*) rekord,lancuch
661 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
662 if (iread.eq.0) return
663 iread=iread+ilen(lancuch)+1
664 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
667 c----------------------------------------------------------------------------
668 subroutine reads(rekord,lancuch,wartosc,default)
670 character*(*) rekord,lancuch,wartosc,default
672 integer ilen,lenlan,lenrec,iread,ireade
678 iread=index(rekord,lancuch(:lenlan)//"=")
679 c print *,"rekord",rekord," lancuch",lancuch
680 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
686 c print *,"iread",iread
687 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
688 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
690 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
692 c print *,"iread",iread
693 if (iread.gt.lenrec) then
698 c print *,"ireade",ireade
699 do while (ireade.lt.lenrec .and.
700 & .not.iblnk(rekord(ireade:ireade)))
703 wartosc=rekord(iread:ireade)
706 c----------------------------------------------------------------------------
707 subroutine multreads(rekord,lancuch,tablica,dim,default)
710 character*(*) rekord,lancuch,tablica(dim),default
712 integer ilen,lenlan,lenrec,iread,ireade
721 iread=index(rekord,lancuch(:lenlan)//"=")
722 c print *,"rekord",rekord," lancuch",lancuch
723 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
724 if (iread.eq.0) return
727 c print *,"iread",iread
728 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
729 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
731 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
733 c print *,"iread",iread
734 if (iread.gt.lenrec) return
736 c print *,"ireade",ireade
737 do while (ireade.lt.lenrec .and.
738 & .not.iblnk(rekord(ireade:ireade)))
741 tablica(i)=rekord(iread:ireade)
745 c----------------------------------------------------------------------------
746 subroutine split_string(rekord,tablica,dim,nsub)
748 integer dim,nsub,i,ii,ll,kk
749 character*(*) tablica(dim)
760 C Find the start of term name
762 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
765 C Parse the name into TABLICA(i) until blank found
766 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
768 tablica(i)(kk:kk)=rekord(ii:ii)
771 if (kk.gt.0) nsub=nsub+1
776 c--------------------------------------------------------------------------------
777 integer function iroof(n,m)
779 if (ii*m .lt. n) ii=ii+1