Fixed the "***" issue in PDB output format by shifting the center of each molecule...
[unres.git] / source / wham / src / readrtns.F
1       subroutine read_general_data(*)
2       implicit none
3       include "DIMENSIONS"
4       include "DIMENSIONS.ZSCOPT"
5       include "DIMENSIONS.FREE"
6       include "COMMON.TORSION"
7       include "COMMON.INTERACT"
8       include "COMMON.IOUNITS"
9       include "COMMON.TIME1"
10       include "COMMON.PROT"
11       include "COMMON.PROTFILES"
12       include "COMMON.CHAIN"
13       include "COMMON.NAMES"
14       include "COMMON.FFIELD"
15       include "COMMON.ENEPS"
16       include "COMMON.WEIGHTS"
17       include "COMMON.FREE"
18       include "COMMON.CONTROL"
19       include "COMMON.ENERGIES"
20       include "COMMON.SBRIDGE"
21       character*800 controlcard
22       integer i,j,k,ii,n_ene_found
23       integer ind,itype1,itype2,itypf,itypsc,itypp
24       integer ilen
25       external ilen
26       character*16 ucase
27       character*16 key
28       external ucase
29
30       call card_concat(controlcard,.true.)
31       call readi(controlcard,"N_ENE",n_ene,max_ene)
32       if (n_ene.gt.max_ene) then
33         write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
34      &    max_ene
35         return1
36       endif
37       call readi(controlcard,"NPARMSET",nparmset,1)
38       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
39       call readi(controlcard,"IPARMPRINT",iparmprint,1)
40       write (iout,*) "PARMPRINT",iparmprint
41       if (nparmset.gt.max_parm) then
42         write (iout,*) "Error: parameter out of range: NPARMSET",
43      &    nparmset, Max_Parm
44         return1
45       endif
46       call readi(controlcard,"MAXIT",maxit,5000)
47       call reada(controlcard,"FIMIN",fimin,1.0d-3)
48       call readi(controlcard,"ENSEMBLES",ensembles,0)
49       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
50       write (iout,*) "Number of energy parameter sets",nparmset
51       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
52       write (iout,*) "MaxSlice",MaxSlice
53       call readi(controlcard,"NSLICE",nslice,1)
54       call flush(iout)
55       if (nslice.gt.MaxSlice) then
56         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
57      &    MaxSlice
58         return1
59       endif
60       write (iout,*) "Frequency of storing conformations",
61      & (isampl(i),i=1,nparmset)
62       write (iout,*) "Maxit",maxit," Fimin",fimin
63       call readi(controlcard,"NQ",nQ,1)
64       if (nQ.gt.MaxQ) then
65         write (iout,*) "Error: parameter out of range: NQ",nq,
66      &    maxq
67         return1
68       endif
69       indpdb=0
70       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
71       call reada(controlcard,"DELTA",delta,1.0d-2)
72       call readi(controlcard,"EINICHECK",einicheck,2)
73       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
74       call readi(controlcard,"NGRIDT",NGridT,400)
75       call reada(controlcard,"STARTGRIDT",StartGridT,200.0d0)
76       call reada(controlcard,"DELTA_T",Delta_T,1.0d0)
77       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
78       call readi(controlcard,"RESCALE",rescale_mode,1)
79       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
80       write (iout,*) "delta",delta
81       write (iout,*) "einicheck",einicheck
82       write (iout,*) "rescale_mode",rescale_mode
83       call flush(iout)
84       bxfile=index(controlcard,"BXFILE").gt.0
85       cxfile=index(controlcard,"CXFILE").gt.0
86       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
87      & bxfile=.true.
88       histfile=index(controlcard,"HISTFILE").gt.0
89       histout=index(controlcard,"HISTOUT").gt.0
90       entfile=index(controlcard,"ENTFILE").gt.0
91       zscfile=index(controlcard,"ZSCFILE").gt.0
92       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
93       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
94       write (iout,*) "with_dihed_constr ",with_dihed_constr,
95      & " CONSTR_DIST",constr_dist
96       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
97       write (iout,*) "with_homology_constr ",with_dihed_constr,
98      & " CONSTR_HOMOLOGY",constr_homology
99       refstr = index(controlcard,'REFSTR').gt.0
100       pdbref = index(controlcard,'PDBREF').gt.0
101       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
102 C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
103 C disulfide bond. Note that in conterary to dynamics this in
104 C CONTROLCARD. The bond are read in molread_zs.F 
105       call flush(iout)
106       return
107       end
108 c------------------------------------------------------------------------------
109       subroutine read_efree(*)
110 C
111 C Read molecular data
112 C
113       implicit none
114       include 'DIMENSIONS'
115       include 'DIMENSIONS.ZSCOPT'
116       include 'DIMENSIONS.COMPAR'
117       include 'DIMENSIONS.FREE'
118       include 'COMMON.IOUNITS'
119       include 'COMMON.TIME1'
120       include 'COMMON.SBRIDGE'
121       include 'COMMON.CONTROL'
122       include 'COMMON.CHAIN'
123       include 'COMMON.HEADER'
124       include 'COMMON.GEO'
125       include 'COMMON.FREE'
126       character*320 controlcard,ucase
127       integer iparm,ib,i,j,npars
128       integer ilen
129       external ilen
130      
131       if (hamil_rep) then
132         npars=1
133       else
134         npars=nParmSet
135       endif
136
137       do iparm=1,npars
138
139       call card_concat(controlcard,.true.)
140       call readi(controlcard,'NT',nT_h(iparm),1)
141       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
142       call flush(iout)
143       if (nT_h(iparm).gt.MaxT_h) then
144         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),
145      &    MaxT_h
146         return1
147       endif
148       replica(iparm)=index(controlcard,"REPLICA").gt.0
149       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
150       if (umbrella(iparm) .and. homol_nset.gt.1) then
151         umbrella(iparm) = .false.
152         write (iout,*) 
153      & "Replica in homology restraints weights UMBRELLA ignored,",iparm
154       endif
155       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
156       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
157      &  replica(iparm)," umbrella ",umbrella(iparm),
158      &  " read_iset",read_iset(iparm)
159       call flush(iout)
160       do ib=1,nT_h(iparm)
161         call card_concat(controlcard,.true.)
162         call readi(controlcard,'NR',nR(ib,iparm),1)
163         if (umbrella(iparm)) then
164           nRR(ib,iparm)=1
165         else
166           nRR(ib,iparm)=nR(ib,iparm)
167         endif
168         if (nR(ib,iparm).gt.MaxR) then
169           write (iout,*)  "Error: parameter out of range: NR",
170      &      nR(ib,iparm),MaxR
171         return1
172         endif
173         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
174         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
175         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
176      &    0.0d0)
177         do i=1,nR(ib,iparm)
178           call card_concat(controlcard,.true.)
179           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
180      &      100.0d0)
181           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
182      &      0.0d0)
183         enddo
184       enddo
185       do ib=1,nT_h(iparm)
186         write (iout,*) "ib",ib," beta_h",
187      &    1.0d0/(0.001987*beta_h(ib,iparm))
188         write (iout,*) "nR",nR(ib,iparm)
189         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
190         do i=1,nR(ib,iparm)
191           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
192      &      "q0",(q0(j,i,ib,iparm),j=1,nQ)
193         enddo
194         call flush(iout)
195       enddo
196
197       write (iout,*) "HOMOL_NSET",homol_nset
198       if (homol_nset.gt.1) then
199         write (iout,*) "HOMOL_NSET: nT_h",nT_h(iparm)
200         do ib=1,nT_h(iparm)
201           nR(ib,iparm)=homol_nset
202           write (iout,*) "iparm",iparm," ib",ib," nR",nR(ib,iparm)
203         enddo
204       endif
205
206       enddo
207
208       if (hamil_rep) then
209
210        do iparm=2,nParmSet
211           nT_h(iparm)=nT_h(1)
212          do ib=1,nT_h(iparm)
213            nR(ib,iparm)=nR(ib,1)
214            if (umbrella(iparm)) then
215              nRR(ib,iparm)=1
216            else
217              nRR(ib,iparm)=nR(ib,1)
218            endif
219            beta_h(ib,iparm)=beta_h(ib,1)
220            do i=1,nR(ib,iparm)
221              f(i,ib,iparm)=f(i,ib,1)
222              do j=1,nQ
223                KH(j,i,ib,iparm)=KH(j,i,ib,1) 
224                Q0(j,i,ib,iparm)=Q0(j,i,ib,1) 
225              enddo
226            enddo
227            replica(iparm)=replica(1)
228            umbrella(iparm)=umbrella(1)
229            read_iset(iparm)=read_iset(1)
230          enddo
231        enddo
232         
233       endif
234
235       return
236       end
237 c-----------------------------------------------------------------------------
238       subroutine read_protein_data(*)
239       implicit none
240       include "DIMENSIONS"
241       include "DIMENSIONS.ZSCOPT"
242       include "DIMENSIONS.FREE"
243 #ifdef MPI
244       include "mpif.h"
245       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
246       include "COMMON.MPI"
247 #endif
248       include "COMMON.CHAIN"
249       include "COMMON.IOUNITS"
250       include "COMMON.PROT"
251       include "COMMON.PROTFILES"
252       include "COMMON.NAMES"
253       include "COMMON.FREE"
254       include "COMMON.OBCINKA"
255       character*64 nazwa
256       character*16000 controlcard
257       integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
258       external ilen,iroof
259       if (hamil_rep) then
260         npars=1
261       else
262         npars=nparmset
263       endif
264
265       do iparm=1,npars
266
267 C Read names of files with conformation data.
268       if (replica(iparm)) then
269         nthr = 1
270       else
271         nthr = nT_h(iparm)
272       endif
273       do ib=1,nthr
274       do ii=1,nRR(ib,iparm)
275       write (iout,*) "Parameter set",iparm," temperature",ib,
276      & " window",ii
277       call flush(iout)
278       call card_concat(controlcard,.true.) 
279       write (iout,*) controlcard(:ilen(controlcard))
280       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
281       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
282       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
283       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
284       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
285      & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
286       call reada(controlcard,"TIME_START",
287      &  time_start_collect(ii,ib,iparm),0.0d0)
288       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
289      &  1.0d10)
290       write (iout,*) "rec_start",rec_start(ii,ib,iparm),
291      & " rec_end",rec_end(ii,ib,iparm)
292       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
293      & " time_end",time_end_collect(ii,ib,iparm)
294       call flush(iout)
295       if (replica(iparm)) then
296         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
297         write (iout,*) "Number of trajectories",totraj(ii,iparm)
298         call flush(iout)
299       endif
300       if (nfile_bin(ii,ib,iparm).lt.2 
301      &    .and. nfile_asc(ii,ib,iparm).eq.0
302      &    .and. nfile_cx(ii,ib,iparm).eq.0) then
303         write (iout,*) "Error - no action specified!"
304         return1
305       endif
306       if (nfile_bin(ii,ib,iparm).gt.0) then
307         call card_concat(controlcard,.false.)
308         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
309      &   maxfile_prot,nfile_bin(ii,ib,iparm))
310 #ifdef DEBUG
311         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
312         write(iout,*) (protfiles(i,1,ii,ib,iparm),
313      &    i=1,nfile_bin(ii,ib,iparm))
314 #endif
315       endif
316       if (nfile_asc(ii,ib,iparm).gt.0) then
317         call card_concat(controlcard,.false.)
318         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
319      &   maxfile_prot,nfile_asc(ii,ib,iparm))
320 #ifdef DEBUG
321         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
322         write(iout,*) (protfiles(i,2,ii,ib,iparm),
323      &    i=1,nfile_asc(ii,ib,iparm))
324 #endif
325       else if (nfile_cx(ii,ib,iparm).gt.0) then
326         call card_concat(controlcard,.false.)
327         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
328      &   maxfile_prot,nfile_cx(ii,ib,iparm))
329 #ifdef DEBUG
330         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
331         write(iout,*) (protfiles(i,2,ii,ib,iparm),
332      &    i=1,nfile_cx(ii,ib,iparm))
333 #endif
334       endif
335       call flush(iout)
336       enddo
337       enddo
338
339       enddo
340
341       return
342       end
343 c-------------------------------------------------------------------------------
344       subroutine opentmp(islice,iunit,bprotfile_temp)
345       implicit none
346       include "DIMENSIONS"
347       include "DIMENSIONS.ZSCOPT"
348       include "DIMENSIONS.FREE"
349 #ifdef MPI
350       include "mpif.h"
351       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
352       include "COMMON.MPI"
353 #endif
354       include "COMMON.IOUNITS"
355       include "COMMON.PROTFILES"
356       include "COMMON.PROT"
357       include "COMMON.FREE"
358       character*64 bprotfile_temp
359       character*3 liczba,liczba2
360       character*2 liczba1
361       integer iunit,islice
362       integer ilen,iroof
363       external ilen,iroof
364       logical lerr
365
366       write (liczba1,'(bz,i2.2)') islice
367       write (liczba,'(bz,i3.3)') me
368 #ifdef MPI
369 c      write (iout,*) "separate_parset ",separate_parset,
370 c     &  " myparm",myparm
371       if (separate_parset) then
372       write (liczba2,'(bz,i3.3)') myparm
373       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
374      &  prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
375       open (iunit,file=bprotfile_temp,status="unknown",
376      &    form="unformatted",access="direct",recl=lenrec)
377       else
378       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
379      &  prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
380       open (iunit,file=bprotfile_temp,status="unknown",
381      &    form="unformatted",access="direct",recl=lenrec)
382       endif
383 #else
384       bprotfile_temp = scratchdir(:ilen(scratchdir))//
385      &  "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
386       open (iunit,file=bprotfile_temp,status="unknown",
387      &    form="unformatted",access="direct",recl=lenrec)
388 #endif      
389 c      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
390 c     &  bprotfile_temp
391 c      call flush(iout)
392       return
393       end
394 c-------------------------------------------------------------------------------
395       subroutine read_database(*)
396       implicit none
397       include "DIMENSIONS"
398       include "DIMENSIONS.ZSCOPT"
399       include "DIMENSIONS.FREE"
400 #ifdef MPI
401       include "mpif.h"
402       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
403       include "COMMON.MPI"
404 #endif
405       include "COMMON.CHAIN"
406       include "COMMON.IOUNITS"
407       include "COMMON.PROTFILES"
408       include "COMMON.NAMES"
409       include "COMMON.VAR"
410       include "COMMON.GEO"
411       include "COMMON.ENEPS"
412       include "COMMON.PROT"
413       include "COMMON.INTERACT"
414       include "COMMON.FREE"
415       include "COMMON.SBRIDGE"
416       include "COMMON.OBCINKA"
417       real*4 csingle(3,maxres2)
418       character*64 nazwa,bprotfile_temp
419       character*3 liczba
420       character*2 liczba1
421       integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
422      &  ll(maxslice),mm(maxslice),if
423       integer nrec,nlines,iscor,iunit,islice
424       double precision energ
425       integer ilen,iroof
426       external ilen,iroof
427       double precision rmsdev,energia(0:max_ene),efree,eini,temp
428       double precision prop(maxQ)
429       integer ntot_all(maxslice,0:maxprocs-1)
430       integer iparm,ib,iib,ir,nprop,nthr,npars
431       double precision etot,time
432       integer ixdrf,iret 
433       logical lerr,linit
434
435       lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
436       lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
437       lenrec=lenrec2+8
438       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
439      &  " lenrec2",lenrec2
440
441       do i=1,nQ
442         prop(i)=0.0d0
443       enddo
444       do islice=1,nslice
445         ll(islice)=0
446         mm(islice)=0
447       enddo
448       write (iout,*) "nparmset",nparmset
449       if (hamil_rep) then
450         npars=1
451       else
452         npars=nparmset
453       endif
454       do iparm=1,npars
455
456       if (replica(iparm)) then
457         nthr = 1
458       else
459         nthr = nT_h(iparm)
460       endif
461
462       do ib=1,nthr
463       do iR=1,nRR(ib,iparm)
464
465       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
466       do islice=1,nslice
467         jj(islice)=0
468         kk(islice)=0
469       enddo
470
471       IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
472 c Read conformations from binary DA files (one per batch) and write them to 
473 c a binary DA scratchfile.
474         write (liczba,'(bz,i3.3)') me
475         do if=1,nfile_bin(iR,ib,iparm)
476           nazwa=protfiles(if,1,iR,ib,iparm)
477      &     (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
478           open (ientin,file=nazwa,status="old",form="unformatted",
479      &     access="direct",recl=lenrec2,err=1111)
480           ii=0
481           do islice=1,nslice
482             call opentmp(islice,ientout,bprotfile_temp)
483             call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
484      &        mm(islice),iR,ib,iparm)
485             close(ientout)
486           enddo
487           close(ientin)
488         enddo
489       ENDIF ! NFILE_BIN>0
490 c
491       IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
492 c Read conformations from multiple ASCII int files and write them to a binary
493 c DA scratchfile.
494         do if=1,nfile_asc(iR,ib,iparm)
495           nazwa=protfiles(if,2,iR,ib,iparm)
496      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
497           open(unit=ientin,file=nazwa,status='old',err=1111)
498           write(iout,*) "reading ",nazwa(:ilen(nazwa))
499           ii=0
500           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
501         enddo ! if
502       ENDIF
503       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
504 c Read conformations from cx files and write them to a binary
505 c DA scratchfile.
506         do if=1,nfile_cx(iR,ib,iparm)
507           nazwa=protfiles(if,2,iR,ib,iparm)
508      &     (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
509           write(iout,*) "reading ",nazwa(:ilen(nazwa))
510           ii=0
511           print *,"Calling cxread"
512           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
513      &       *1111)
514           close(ientout)
515           write (iout,*) "exit cxread"
516           call flush(iout)
517         enddo
518       ENDIF
519
520       do islice=1,nslice
521         stot(islice)=stot(islice)+jj(islice)
522       enddo
523
524       enddo
525       enddo
526       write (iout,*) "IPARM",iparm
527       enddo
528
529       if (nslice.eq.1) then
530 #ifdef MPI
531         write (liczba,'(bz,i3.3)') me
532         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
533      &    prefix(:ilen(prefix))//liczba//".xbin.tmp"
534 #else
535         bprotfile_temp = scratchdir(:ilen(scratchdir))//
536      &     "/"//prefix(:ilen(prefix))//".xbin.tmp"
537 #endif
538         write(iout,*) mm(1)," conformations read",ll(1),
539      &    " conformations written to ", 
540      &    bprotfile_temp(:ilen(bprotfile_temp))
541       else
542         do islice=1,nslice
543           write (liczba1,'(bz,i2.2)') islice
544 #ifdef MPI
545           write (liczba,'(bz,i3.3)') me
546           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
547      &      prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
548 #else
549           bprotfile_temp = scratchdir(:ilen(scratchdir))//
550      &       "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
551 #endif
552           write(iout,*) mm(islice)," conformations read",ll(islice),
553      &    " conformations written to ", 
554      &    bprotfile_temp(:ilen(bprotfile_temp))
555         enddo
556       endif
557
558 #ifdef MPI
559 c Check if everyone has the same number of conformations
560       call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
561      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
562       lerr=.false.
563       do i=0,nprocs-1
564         if (i.ne.me) then
565           do islice=1,nslice
566           if (stot(islice).ne.ntot_all(islice,i)) then
567             write (iout,*) "Number of conformations at processor",i,
568      &       " differs from that at processor",me,
569      &       stot(islice),ntot_all(islice,i)," slice",islice
570             lerr = .true.
571           endif
572           enddo
573         endif
574       enddo 
575       if (lerr) then
576         write (iout,*)
577         write (iout,*) "Numbers of conformations read by processors"
578         write (iout,*)
579         do i=0,nprocs-1
580           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
581         enddo
582         write (iout,*) "Calculation terminated."
583         call flush(iout)
584         return1
585       endif
586       do islice=1,nslice
587         ntot(islice)=stot(islice)
588       enddo
589       return
590 #endif
591  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
592       call flush(iout)
593       return1
594       end
595 c------------------------------------------------------------------------------
596       subroutine card_concat(card,to_upper)
597       implicit none
598       include 'DIMENSIONS.ZSCOPT'
599       include "COMMON.IOUNITS"
600       character*(*) card
601       character*80 karta,ucase
602       logical to_upper
603       integer ilen
604       external ilen
605       read (inp,'(a)') karta
606       if (to_upper) karta=ucase(karta)
607       card=' '
608       do while (karta(80:80).eq.'&')
609         card=card(:ilen(card)+1)//karta(:79)
610         read (inp,'(a)') karta
611         if (to_upper) karta=ucase(karta)
612       enddo
613       card=card(:ilen(card)+1)//karta
614       return
615       end
616 c------------------------------------------------------------------------------
617       subroutine readi(rekord,lancuch,wartosc,default)
618       implicit none
619       character*(*) rekord,lancuch
620       integer wartosc,default
621       integer ilen,iread
622       external ilen
623       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
624       if (iread.eq.0) then
625         wartosc=default
626         return
627       endif
628       iread=iread+ilen(lancuch)+1
629       read (rekord(iread:),*) wartosc
630       return
631       end
632 c----------------------------------------------------------------------------
633       subroutine reada(rekord,lancuch,wartosc,default)
634       implicit none
635       character*(*) rekord,lancuch
636       character*80 aux
637       double precision wartosc,default
638       integer ilen,iread
639       external ilen
640       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
641       if (iread.eq.0) then
642         wartosc=default
643         return
644       endif
645       iread=iread+ilen(lancuch)+1
646       read (rekord(iread:),*) wartosc
647       return
648       end
649 c----------------------------------------------------------------------------
650       subroutine multreadi(rekord,lancuch,tablica,dim,default)
651       implicit none
652       integer dim,i
653       integer tablica(dim),default
654       character*(*) rekord,lancuch
655       character*80 aux
656       integer ilen,iread
657       external ilen
658       do i=1,dim
659         tablica(i)=default
660       enddo
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)
665    10 return
666       end
667 c----------------------------------------------------------------------------
668       subroutine multreada(rekord,lancuch,tablica,dim,default)
669       implicit none
670       integer dim,i
671       double precision tablica(dim),default
672       character*(*) rekord,lancuch
673       character*80 aux
674       integer ilen,iread
675       external ilen
676       do i=1,dim
677         tablica(i)=default
678       enddo
679       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
680       if (iread.eq.0) return
681       iread=iread+ilen(lancuch)+1
682       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
683    10 return
684       end
685 c----------------------------------------------------------------------------
686       subroutine reads(rekord,lancuch,wartosc,default)
687       implicit none
688       character*(*) rekord,lancuch,wartosc,default
689       character*80 aux
690       integer ilen,lenlan,lenrec,iread,ireade
691       external ilen
692       logical iblnk
693       external iblnk
694       lenlan=ilen(lancuch)
695       lenrec=ilen(rekord)
696       iread=index(rekord,lancuch(:lenlan)//"=")
697 c      print *,"rekord",rekord," lancuch",lancuch
698 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
699       if (iread.eq.0) then
700         wartosc=default
701         return
702       endif
703       iread=iread+lenlan+1
704 c      print *,"iread",iread
705 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
706       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
707         iread=iread+1
708 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
709       enddo
710 c      print *,"iread",iread
711       if (iread.gt.lenrec) then
712          wartosc=default
713         return
714       endif
715       ireade=iread+1
716 c      print *,"ireade",ireade
717       do while (ireade.lt.lenrec .and.
718      &   .not.iblnk(rekord(ireade:ireade)))
719         ireade=ireade+1
720       enddo
721       wartosc=rekord(iread:ireade)
722       return
723       end
724 c----------------------------------------------------------------------------
725       subroutine multreads(rekord,lancuch,tablica,dim,default)
726       implicit none
727       integer dim,i
728       character*(*) rekord,lancuch,tablica(dim),default
729       character*80 aux
730       integer ilen,lenlan,lenrec,iread,ireade
731       external ilen
732       logical iblnk
733       external iblnk
734       do i=1,dim
735         tablica(i)=default
736       enddo
737       lenlan=ilen(lancuch)
738       lenrec=ilen(rekord)
739       iread=index(rekord,lancuch(:lenlan)//"=")
740 c      print *,"rekord",rekord," lancuch",lancuch
741 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
742       if (iread.eq.0) return
743       iread=iread+lenlan+1
744       do i=1,dim
745 c      print *,"iread",iread
746 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
747       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
748         iread=iread+1
749 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
750       enddo
751 c      print *,"iread",iread
752       if (iread.gt.lenrec) return
753       ireade=iread+1
754 c      print *,"ireade",ireade
755       do while (ireade.lt.lenrec .and.
756      &   .not.iblnk(rekord(ireade:ireade)))
757         ireade=ireade+1
758       enddo
759       tablica(i)=rekord(iread:ireade)
760       iread=ireade+1
761       enddo
762       end
763 c----------------------------------------------------------------------------
764       subroutine split_string(rekord,tablica,dim,nsub)
765       implicit none
766       integer dim,nsub,i,ii,ll,kk
767       character*(*) tablica(dim)
768       character*(*) rekord
769       integer ilen
770       external ilen
771       do i=1,dim
772         tablica(i)=" "
773       enddo
774       ii=1
775       ll = ilen(rekord)
776       nsub=0
777       do i=1,dim
778 C Find the start of term name
779         kk = 0
780         do while (ii.le.ll .and. rekord(ii:ii).eq." ") 
781           ii = ii+1
782         enddo
783 C Parse the name into TABLICA(i) until blank found
784         do while (ii.le.ll .and. rekord(ii:ii).ne." ") 
785           kk = kk+1
786           tablica(i)(kk:kk)=rekord(ii:ii)
787           ii = ii+1
788         enddo
789         if (kk.gt.0) nsub=nsub+1
790         if (ii.gt.ll) return
791       enddo
792       return
793       end
794 c--------------------------------------------------------------------------------
795       integer function iroof(n,m)
796       ii = n/m
797       if (ii*m .lt. n) ii=ii+1
798       iroof = ii
799       return
800       end