9d444cfaaa0b9eaec0058dafa716f32115fed104
[unres.git] / readrtns_min.F
1       subroutine readrtns
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4 #ifdef MPI
5       include 'mpif.h'
6 #endif
7       include 'COMMON.SETUP'
8       include 'COMMON.CONTROL'
9       include 'COMMON.SBRIDGE'
10       include 'COMMON.IOUNITS'
11       logical file_exist
12 C Read force-field parameters except weights
13       call parmread
14 C Read job setup parameters
15       call read_control
16 C Read control parameters for energy minimzation if required
17       if (minim) call read_minim
18 C Read MCM control parameters if required
19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 c      if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 c      if (modecalc.eq.14) then 
24 c         call read_MDpar
25 c         call read_REMDpar
26 c      endif
27 C Read MUCA control parameters if required
28 c      if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 c      if (modecalc.eq.8) then
32 c       inquire (file="fort.40",exist=file_exist)
33 c       if (.not.file_exist) call csaread
34 c      endif 
35 cfmc      if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
38       call molread
39 C Print restraint information
40 #ifdef MPI
41       if (.not. out1file .or. me.eq.king) then
42 #endif
43       if (nhpb.gt.nss) 
44      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45      & " distance constraints have been imposed"
46       do i=nss+1,nhpb
47         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
48       enddo
49 #ifdef MPI
50       endif
51 #endif
52 c      print *,"Processor",myrank," leaves READRTNS"
53       return
54       end
55 C-------------------------------------------------------------------------------
56       subroutine read_control
57 C
58 C Read contorl data
59 C
60       implicit real*8 (a-h,o-z)
61       include 'DIMENSIONS'
62 #ifdef MP
63       include 'mpif.h'
64       logical OKRandom, prng_restart
65       real*8  r1
66 #endif
67       include 'COMMON.IOUNITS'
68       include 'COMMON.TIME1'
69       include 'COMMON.SBRIDGE'
70       include 'COMMON.CONTROL'
71       include 'COMMON.MCM'
72       include 'COMMON.HEADER'
73       include 'COMMON.CHAIN'
74       include 'COMMON.FFIELD'
75       include 'COMMON.SETUP'
76       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
77       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
78       character*80 ucase
79       character*320 controlcard
80
81       nglob_csa=0
82       eglob_csa=1d99
83       nmin_csa=0
84       read (INP,'(a)') titel
85       call card_concat(controlcard)
86 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
87 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
88       call reada(controlcard,'SEED',seed,0.0D0)
89       call random_init(seed)
90 C Set up the time limit (caution! The time must be input in minutes!)
91       read_cart=index(controlcard,'READ_CART').gt.0
92       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
93       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
94       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
95       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
96       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
97       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
98       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
99       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
100       call reada(controlcard,'DRMS',drms,0.1D0)
101       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
102        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc 
103        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 
104        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max 
105        write (iout,'(a,f10.1)')'DRMS    = ',drms 
106        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm 
107        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
108       endif
109       call readi(controlcard,'NZ_START',nz_start,0)
110       call readi(controlcard,'NZ_END',nz_end,0)
111       call readi(controlcard,'IZ_SC',iz_sc,0)
112       timlim=60.0D0*timlim
113       safety = 60.0d0*safety
114       timem=timlim
115       modecalc=0
116       call reada(controlcard,"T_BATH",t_bath,300.0d0)
117       minim=(index(controlcard,'MINIMIZE').gt.0)
118       dccart=(index(controlcard,'CART').gt.0)
119       overlapsc=(index(controlcard,'OVERLAP').gt.0)
120       overlapsc=.not.overlapsc
121       searchsc=(index(controlcard,'SEARCHSC').gt.0)
122       sideadd=(index(controlcard,'SIDEADD').gt.0)
123       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
124       outpdb=(index(controlcard,'PDBOUT').gt.0)
125       outmol2=(index(controlcard,'MOL2OUT').gt.0)
126       pdbref=(index(controlcard,'PDBREF').gt.0)
127       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
128       indpdb=index(controlcard,'PDBSTART')
129       extconf=(index(controlcard,'EXTCONF').gt.0)
130       call readi(controlcard,'IPRINT',iprint,0)
131       call readi(controlcard,'MAXGEN',maxgen,10000)
132       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
133       call readi(controlcard,"KDIAG",kdiag,0)
134       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
135       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
136      & write (iout,*) "RESCALE_MODE",rescale_mode
137       split_ene=index(controlcard,'SPLIT_ENE').gt.0
138       if (index(controlcard,'REGULAR').gt.0.0D0) then
139         call reada(controlcard,'WEIDIS',weidis,0.1D0)
140         modecalc=1
141         refstr=.true.
142       endif
143       if (index(controlcard,'CHECKGRAD').gt.0) then
144         modecalc=5
145         if (index(controlcard,'CART').gt.0) then
146           icheckgrad=1
147         elseif (index(controlcard,'CARINT').gt.0) then
148           icheckgrad=2
149         else
150           icheckgrad=3
151         endif
152       elseif (index(controlcard,'THREAD').gt.0) then
153         modecalc=2
154         call readi(controlcard,'THREAD',nthread,0)
155         if (nthread.gt.0) then
156           call reada(controlcard,'WEIDIS',weidis,0.1D0)
157         else
158           if (fg_rank.eq.0)
159      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
160           stop 'Error termination in Read_Control.'
161         endif
162       else if (index(controlcard,'MCMA').gt.0) then
163         modecalc=3
164       else if (index(controlcard,'MCEE').gt.0) then
165         modecalc=6
166       else if (index(controlcard,'MULTCONF').gt.0) then
167         modecalc=4
168       else if (index(controlcard,'MAP').gt.0) then
169         modecalc=7
170         call readi(controlcard,'MAP',nmap,0)
171       else if (index(controlcard,'CSA').gt.0) then
172         modecalc=8
173 crc      else if (index(controlcard,'ZSCORE').gt.0) then
174 crc   
175 crc  ZSCORE is rm from UNRES, modecalc=9 is available
176 crc
177 crc        modecalc=9
178 cfcm      else if (index(controlcard,'MCMF').gt.0) then
179 cfmc        modecalc=10
180       else if (index(controlcard,'SOFTREG').gt.0) then
181         modecalc=11
182       else if (index(controlcard,'CHECK_BOND').gt.0) then
183         modecalc=-1
184       else if (index(controlcard,'TEST').gt.0) then
185         modecalc=-2
186       else if (index(controlcard,'MD').gt.0) then
187         modecalc=12
188       else if (index(controlcard,'RE ').gt.0) then
189         modecalc=14
190       endif
191
192       lmuca=index(controlcard,'MUCA').gt.0
193       call readi(controlcard,'MUCADYN',mucadyn,0)      
194       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
195       if (lmuca .and. (me.eq.king .or. .not.out1file )) 
196      & then
197        write (iout,*) 'MUCADYN=',mucadyn
198        write (iout,*) 'MUCASMOOTH=',muca_smooth
199       endif
200
201       iscode=index(controlcard,'ONE_LETTER')
202       indphi=index(controlcard,'PHI')
203       indback=index(controlcard,'BACK')
204       iranconf=index(controlcard,'RAND_CONF')
205       i2ndstr=index(controlcard,'USE_SEC_PRED')
206       gradout=index(controlcard,'GRADOUT').gt.0
207       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
208       
209       if(me.eq.king.or..not.out1file)
210      & write (iout,'(2a)') diagmeth(kdiag),
211      &  ' routine used to diagonalize matrices.'
212       return
213       end
214 c------------------------------------------------------------------------------
215       subroutine molread
216 C
217 C Read molecular data.
218 C
219       implicit real*8 (a-h,o-z)
220       include 'DIMENSIONS'
221 #ifdef MPI
222       include 'mpif.h'
223       integer error_msg
224 #endif
225       include 'COMMON.IOUNITS'
226       include 'COMMON.GEO'
227       include 'COMMON.VAR'
228       include 'COMMON.INTERACT'
229       include 'COMMON.LOCAL'
230       include 'COMMON.NAMES'
231       include 'COMMON.CHAIN'
232       include 'COMMON.FFIELD'
233       include 'COMMON.SBRIDGE'
234       include 'COMMON.HEADER'
235       include 'COMMON.CONTROL'
236 c      include 'COMMON.DBASE'
237 c      include 'COMMON.THREAD'
238       include 'COMMON.CONTACTS'
239       include 'COMMON.TORCNSTR'
240       include 'COMMON.TIME1'
241       include 'COMMON.BOUNDS'
242 c      include 'COMMON.MD'
243 c      include 'COMMON.REMD'
244       include 'COMMON.SETUP'
245       character*4 sequence(maxres)
246       integer rescode
247       double precision x(maxvar)
248       character*256 pdbfile
249       character*320 weightcard
250       character*80 weightcard_t,ucase
251       dimension itype_pdb(maxres)
252       common /pizda/ itype_pdb
253       logical seq_comp,fail
254       double precision energia(0:n_ene)
255       integer ilen
256       external ilen
257 C
258 C Body
259 C
260 C Read weights of the subsequent energy terms.
261
262        call card_concat(weightcard)
263        call reada(weightcard,'WLONG',wlong,1.0D0)
264        call reada(weightcard,'WSC',wsc,wlong)
265        call reada(weightcard,'WSCP',wscp,wlong)
266        call reada(weightcard,'WELEC',welec,1.0D0)
267        call reada(weightcard,'WVDWPP',wvdwpp,welec)
268        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
269        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
270        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
271        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
272        call reada(weightcard,'WTURN3',wturn3,1.0D0)
273        call reada(weightcard,'WTURN4',wturn4,1.0D0)
274        call reada(weightcard,'WTURN6',wturn6,1.0D0)
275        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
276        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
277        call reada(weightcard,'WBOND',wbond,1.0D0)
278        call reada(weightcard,'WTOR',wtor,1.0D0)
279        call reada(weightcard,'WTORD',wtor_d,1.0D0)
280        call reada(weightcard,'WANG',wang,1.0D0)
281        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
282        call reada(weightcard,'SCAL14',scal14,0.4D0)
283        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
284        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
285        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
286        call reada(weightcard,'TEMP0',temp0,300.0d0)
287        if (index(weightcard,'SOFT').gt.0) ipot=6
288 C 12/1/95 Added weight for the multi-body term WCORR
289        call reada(weightcard,'WCORRH',wcorr,1.0D0)
290        if (wcorr4.gt.0.0d0) wcorr=wcorr4
291        weights(1)=wsc
292        weights(2)=wscp
293        weights(3)=welec
294        weights(4)=wcorr
295        weights(5)=wcorr5
296        weights(6)=wcorr6
297        weights(7)=wel_loc
298        weights(8)=wturn3
299        weights(9)=wturn4
300        weights(10)=wturn6
301        weights(11)=wang
302        weights(12)=wscloc
303        weights(13)=wtor
304        weights(14)=wtor_d
305        weights(15)=wstrain
306        weights(16)=wvdwpp
307        weights(17)=wbond
308        weights(18)=scal14
309        weights(21)=wsccor
310
311       if(me.eq.king.or..not.out1file)
312      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
313      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
314      &  wturn4,wturn6
315    10 format (/'Energy-term weights (unscaled):'//
316      & 'WSCC=   ',f10.6,' (SC-SC)'/
317      & 'WSCP=   ',f10.6,' (SC-p)'/
318      & 'WELEC=  ',f10.6,' (p-p electr)'/
319      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
320      & 'WBOND=  ',f10.6,' (stretching)'/
321      & 'WANG=   ',f10.6,' (bending)'/
322      & 'WSCLOC= ',f10.6,' (SC local)'/
323      & 'WTOR=   ',f10.6,' (torsional)'/
324      & 'WTORD=  ',f10.6,' (double torsional)'/
325      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
326      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
327      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
328      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
329      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
330      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
331      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
332      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
333      & 'WTURN6= ',f10.6,' (turns, 6th order)')
334       if(me.eq.king.or..not.out1file)then
335        if (wcorr4.gt.0.0d0) then
336         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
337      &   'between contact pairs of peptide groups'
338         write (iout,'(2(a,f5.3/))') 
339      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
340      &  'Range of quenching the correlation terms:',2*delt_corr 
341        else if (wcorr.gt.0.0d0) then
342         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
343      &   'between contact pairs of peptide groups'
344        endif
345        write (iout,'(a,f8.3)') 
346      &  'Scaling factor of 1,4 SC-p interactions:',scal14
347        write (iout,'(a,f8.3)') 
348      &  'General scaling factor of SC-p interactions:',scalscp
349       endif
350       r0_corr=cutoff_corr-delt_corr
351       do i=1,20
352         aad(i,1)=scalscp*aad(i,1)
353         aad(i,2)=scalscp*aad(i,2)
354         bad(i,1)=scalscp*bad(i,1)
355         bad(i,2)=scalscp*bad(i,2)
356       enddo
357 c      call rescale_weights(t_bath)
358       if(me.eq.king.or..not.out1file)
359      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
360      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
361      &  wturn4,wturn6
362    22 format (/'Energy-term weights (scaled):'//
363      & 'WSCC=   ',f10.6,' (SC-SC)'/
364      & 'WSCP=   ',f10.6,' (SC-p)'/
365      & 'WELEC=  ',f10.6,' (p-p electr)'/
366      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
367      & 'WBOND=  ',f10.6,' (stretching)'/
368      & 'WANG=   ',f10.6,' (bending)'/
369      & 'WSCLOC= ',f10.6,' (SC local)'/
370      & 'WTOR=   ',f10.6,' (torsional)'/
371      & 'WTORD=  ',f10.6,' (double torsional)'/
372      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
373      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
374      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
375      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
376      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
377      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
378      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
379      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
380      & 'WTURN6= ',f10.6,' (turns, 6th order)')
381       if(me.eq.king.or..not.out1file)
382      & write (iout,*) "Reference temperature for weights calculation:",
383      &  temp0
384       call reada(weightcard,"D0CM",d0cm,3.78d0)
385       call reada(weightcard,"AKCM",akcm,15.1d0)
386       call reada(weightcard,"AKTH",akth,11.0d0)
387       call reada(weightcard,"AKCT",akct,12.0d0)
388       call reada(weightcard,"V1SS",v1ss,-1.08d0)
389       call reada(weightcard,"V2SS",v2ss,7.61d0)
390       call reada(weightcard,"V3SS",v3ss,13.7d0)
391       call reada(weightcard,"EBR",ebr,-5.50D0)
392       if(me.eq.king.or..not.out1file) then
393        write (iout,*) "Parameters of the SS-bond potential:"
394        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
395      & " AKCT",akct
396        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
397        write (iout,*) "EBR",ebr
398        print *,'indpdb=',indpdb,' pdbref=',pdbref
399       endif
400       if (indpdb.gt.0 .or. pdbref) then
401         read(inp,'(a)') pdbfile
402         if(me.eq.king.or..not.out1file)
403      &   write (iout,'(2a)') 'PDB data will be read from file ',
404      &   pdbfile(:ilen(pdbfile))
405         open(ipdbin,file=pdbfile,status='old',err=33)
406         goto 34 
407   33    write (iout,'(a)') 'Error opening PDB file.'
408         stop
409   34    continue
410 c        print *,'Begin reading pdb data'
411         call readpdb
412 c        print *,'Finished reading pdb data'
413         if(me.eq.king.or..not.out1file)
414      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
415      &   ' nstart_sup=',nstart_sup
416         do i=1,nres
417           itype_pdb(i)=itype(i)
418         enddo
419         close (ipdbin)
420         nnt=nstart_sup
421         nct=nstart_sup+nsup-1
422 c        call contact(.false.,ncont_ref,icont_ref,co)
423
424         if (sideadd) then 
425 C Following 2 lines for diagnostics; comment out if not needed
426          write (iout,*) "Before sideadd"
427          call intout
428          if(me.eq.king.or..not.out1file)
429      &    write(iout,*)'Adding sidechains'
430          maxsi=1000
431          do i=2,nres-1
432           iti=itype(i)
433           if (iti.ne.10) then
434             nsi=0
435             fail=.true.
436             do while (fail.and.nsi.le.maxsi)
437               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
438               nsi=nsi+1
439             enddo
440             if(fail) write(iout,*)'Adding sidechain failed for res ',
441      &              i,' after ',nsi,' trials'
442           endif
443          enddo
444 C 9/29/12 Adam: Recalculate coordinates with new side chain positions
445          call chainbuild
446         endif  
447 C Following 2 lines for diagnostics; comment out if not needed
448         write (iout,*) "After sideadd"
449         call intout
450       endif
451       if (indpdb.eq.0) then
452 C Read sequence if not taken from the pdb file.
453         read (inp,*) nres
454 c        print *,'nres=',nres
455         if (iscode.gt.0) then
456           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
457         else
458           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
459         endif
460 C Convert sequence to numeric code
461         do i=1,nres
462           itype(i)=rescode(i,sequence(i),iscode)
463         enddo
464 C Assign initial virtual bond lengths
465         do i=2,nres
466           vbld(i)=vbl
467           vbld_inv(i)=vblinv
468         enddo
469         do i=2,nres-1
470           vbld(i+nres)=dsc(itype(i))
471           vbld_inv(i+nres)=dsc_inv(itype(i))
472 c          write (iout,*) "i",i," itype",itype(i),
473 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
474         enddo
475       endif 
476 c      print *,nres
477 c      print '(20i4)',(itype(i),i=1,nres)
478       do i=1,nres
479 #ifdef PROCOR
480         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
481 #else
482         if (itype(i).eq.21) then
483 #endif
484           itel(i)=0
485 #ifdef PROCOR
486         else if (itype(i+1).ne.20) then
487 #else
488         else if (itype(i).ne.20) then
489 #endif
490           itel(i)=1
491         else
492           itel(i)=2
493         endif  
494       enddo
495       if(me.eq.king.or..not.out1file)then
496        write (iout,*) "ITEL"
497        do i=1,nres-1
498          write (iout,*) i,itype(i),itel(i)
499        enddo
500        print *,'Call Read_Bridge.'
501       endif
502       call read_bridge
503 C 8/13/98 Set limits to generating the dihedral angles
504       do i=1,nres
505         phibound(1,i)=-pi
506         phibound(2,i)=pi
507       enddo
508       read (inp,*) ndih_constr
509       if (ndih_constr.gt.0) then
510         read (inp,*) ftors
511         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
512         if(me.eq.king.or..not.out1file)then
513          write (iout,*) 
514      &   'There are',ndih_constr,' constraints on phi angles.'
515          do i=1,ndih_constr
516           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
517          enddo
518         endif
519         do i=1,ndih_constr
520           phi0(i)=deg2rad*phi0(i)
521           drange(i)=deg2rad*drange(i)
522         enddo
523         if(me.eq.king.or..not.out1file)
524      &   write (iout,*) 'FTORS',ftors
525         do i=1,ndih_constr
526           ii = idih_constr(i)
527           phibound(1,ii) = phi0(i)-drange(i)
528           phibound(2,ii) = phi0(i)+drange(i)
529         enddo 
530       endif
531       nnt=1
532 #ifdef MPI
533       if (me.eq.king) then
534 #endif
535        write (iout,'(a)') 'Boundaries in phi angle sampling:'
536        do i=1,nres
537          write (iout,'(a3,i5,2f10.1)') 
538      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
539        enddo
540 #ifdef MP
541       endif
542 #endif
543       nct=nres
544 cd      print *,'NNT=',NNT,' NCT=',NCT
545       if (itype(1).eq.21) nnt=2
546       if (itype(nres).eq.21) nct=nct-1
547       if (pdbref) then
548         if(me.eq.king.or..not.out1file)
549      &   write (iout,'(a,i3)') 'nsup=',nsup
550         nstart_seq=nnt
551         if (nsup.le.(nct-nnt+1)) then
552           do i=0,nct-nnt+1-nsup
553             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
554               nstart_seq=nnt+i
555               goto 111
556             endif
557           enddo
558           write (iout,'(a)') 
559      &            'Error - sequences to be superposed do not match.'
560           stop
561         else
562           do i=0,nsup-(nct-nnt+1)
563             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
564      &      then
565               nstart_sup=nstart_sup+i
566               nsup=nct-nnt+1
567               goto 111
568             endif
569           enddo 
570           write (iout,'(a)') 
571      &            'Error - sequences to be superposed do not match.'
572         endif
573   111   continue
574         if (nsup.eq.0) nsup=nct-nnt
575         if (nstart_sup.eq.0) nstart_sup=nnt
576         if (nstart_seq.eq.0) nstart_seq=nnt
577         if(me.eq.king.or..not.out1file)  
578      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
579      &                 ' nstart_seq=',nstart_seq
580       endif
581 c--- Zscore rms -------
582       if (nz_start.eq.0) nz_start=nnt
583       if (nz_end.eq.0 .and. nsup.gt.0) then
584         nz_end=nnt+nsup-1
585       else if (nz_end.eq.0) then
586         nz_end=nct
587       endif
588       if(me.eq.king.or..not.out1file)then
589        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
590        write (iout,*) 'IZ_SC=',iz_sc
591       endif
592 c----------------------
593       call init_int_table
594       if (refstr) then
595         if (.not.pdbref) then
596           call read_angles(inp,*38)
597           goto 39
598    38     write (iout,'(a)') 'Error reading reference structure.'
599 #ifdef MPI
600           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
601           stop 'Error reading reference structure'
602 #endif
603    39     call chainbuild
604           call setup_var
605 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
606           nstart_sup=nnt
607           nstart_seq=nnt
608           nsup=nct-nnt+1
609           do i=1,2*nres
610             do j=1,3
611               cref(j,i)=c(j,i)
612             enddo
613           enddo
614 c          call contact(.true.,ncont_ref,icont_ref,co)
615         endif
616 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
617         call flush(iout)
618         if (constr_dist.gt.0) call read_dist_constr
619 c        write (iout,*) "After read_dist_constr nhpb",nhpb
620         call hpb_partition
621         if(me.eq.king.or..not.out1file)
622      &   write (iout,*) 'Contact order:',co
623         if (pdbref) then
624         if(me.eq.king.or..not.out1file)
625      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
626         do i=1,ncont_ref
627           do j=1,2
628             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
629           enddo
630           if(me.eq.king.or..not.out1file)
631      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
632      &     icont_ref(1,i),' ',
633      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
634         enddo
635         endif
636       endif
637       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
638      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
639      &    modecalc.ne.10) then
640 C If input structure hasn't been supplied from the PDB file read or generate
641 C initial geometry.
642         if (iranconf.eq.0 .and. .not. extconf) then
643           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
644      &     write (iout,'(a)') 'Initial geometry will be read in.'
645           if (read_cart) then
646             read(inp,'(8f10.5)',end=36,err=36)
647      &       ((c(l,k),l=1,3),k=1,nres),
648      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
649             call int_from_cart1(.false.)
650             do i=1,nres-1
651               do j=1,3
652                 dc(j,i)=c(j,i+1)-c(j,i)
653                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
654               enddo
655             enddo
656             do i=nnt,nct
657               if (itype(i).ne.10) then
658                 do j=1,3
659                   dc(j,i+nres)=c(j,i+nres)-c(j,i) 
660                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
661                 enddo
662               endif
663             enddo
664             return
665           else
666             call read_angles(inp,*36)
667           endif
668           goto 37
669    36     write (iout,'(a)') 'Error reading angle file.'
670 #ifdef MPI
671           call mpi_finalize( MPI_COMM_WORLD,IERR )
672 #endif
673           stop 'Error reading angle file.'
674    37     continue 
675         else if (extconf) then
676          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
677      &    write (iout,'(a)') 'Extended chain initial geometry.'
678          do i=3,nres
679           theta(i)=90d0*deg2rad
680          enddo
681          do i=4,nres
682           phi(i)=180d0*deg2rad
683          enddo
684          do i=2,nres-1
685           alph(i)=110d0*deg2rad
686          enddo
687          do i=2,nres-1
688           omeg(i)=-120d0*deg2rad
689          enddo
690         else
691           if(me.eq.king.or..not.out1file)
692      &     write (iout,'(a)') 'Random-generated initial geometry.'
693
694
695 #ifdef MPI
696           if (me.eq.king  .or. fg_rank.eq.0 .and. (
697      &           modecalc.eq.12 .or. modecalc.eq.14) ) then  
698 #endif
699             do itrial=1,100
700               itmp=1
701               call gen_rand_conf(itmp,*30)
702               goto 40
703    30         write (iout,*) 'Failed to generate random conformation',
704      &          ', itrial=',itrial
705               write (*,*) 'Processor:',me,
706      &          ' Failed to generate random conformation',
707      &          ' itrial=',itrial
708               call intout
709
710 #ifdef AIX
711               call flush_(iout)
712 #else
713               call flush(iout)
714 #endif
715             enddo
716             write (iout,'(a,i3,a)') 'Processor:',me,
717      &        ' error in generating random conformation.'
718             write (*,'(a,i3,a)') 'Processor:',me,
719      &        ' error in generating random conformation.'
720             call flush(iout)
721 #ifdef MPI
722             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)            
723    40       continue
724           endif
725 #else
726           do itrial=1,100
727             itmp=1
728             call gen_rand_conf(itmp,*31)
729             goto 40
730    31       write (iout,*) 'Failed to generate random conformation',
731      &        ', itrial=',itrial
732             write (*,*) 'Failed to generate random conformation',
733      &        ', itrial=',itrial
734           enddo
735           write (iout,'(a,i3,a)') 'Processor:',me,
736      &      ' error in generating random conformation.'
737           write (*,'(a,i3,a)') 'Processor:',me,
738      &      ' error in generating random conformation.'
739           stop
740    40     continue
741 #endif
742         endif
743       elseif (modecalc.eq.4) then
744         read (inp,'(a)') intinname
745         open (intin,file=intinname,status='old',err=333)
746         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
747      &  write (iout,'(a)') 'intinname',intinname
748         write (*,'(a)') 'Processor',myrank,' intinname',intinname
749         goto 334
750   333   write (iout,'(2a)') 'Error opening angle file ',intinname
751 #ifdef MPI 
752         call MPI_Finalize(MPI_COMM_WORLD,IERR)
753 #endif   
754         stop 'Error opening angle file.' 
755   334   continue
756
757       endif 
758 C Generate distance constraints, if the PDB structure is to be regularized. 
759 c      if (nthread.gt.0) then
760 c        call read_threadbase
761 c      endif
762       call setup_var
763       if (me.eq.king .or. .not. out1file)
764      & call intout
765       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
766         write (iout,'(/a,i3,a)') 
767      &  'The chain contains',ns,' disulfide-bridging cysteines.'
768         write (iout,'(20i4)') (iss(i),i=1,ns)
769         write (iout,'(/a/)') 'Pre-formed links are:' 
770         do i=1,nss
771           i1=ihpb(i)-nres
772           i2=jhpb(i)-nres
773           it1=itype(i1)
774           it2=itype(i2)
775           if (me.eq.king.or..not.out1file)
776      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
777      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
778      &    ebr,forcon(i)
779         enddo
780         write (iout,'(a)')
781       endif
782 c       if (i2ndstr.gt.0) call secstrp2dihc
783 c      call geom_to_var(nvar,x)
784 c      call etotal(energia(0))
785 c      call enerprint(energia(0))
786 c      call briefout(0,etot)
787 c      stop
788 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
789 cd    write (iout,'(a)') 'Variable list:'
790 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
791 #ifdef MPI
792       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
793      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
794      &  'Processor',myrank,': end reading molecular data.'
795 #endif
796       return
797       end
798 c--------------------------------------------------------------------------
799       logical function seq_comp(itypea,itypeb,length)
800       implicit none
801       integer length,itypea(length),itypeb(length)
802       integer i
803       do i=1,length
804         if (itypea(i).ne.itypeb(i)) then
805           seq_comp=.false.
806           return
807         endif
808       enddo
809       seq_comp=.true.
810       return
811       end
812 c-----------------------------------------------------------------------------
813       subroutine read_bridge
814 C Read information about disulfide bridges.
815       implicit real*8 (a-h,o-z)
816       include 'DIMENSIONS'
817 #ifdef MPI
818       include 'mpif.h'
819 #endif
820       include 'COMMON.IOUNITS'
821       include 'COMMON.GEO'
822       include 'COMMON.VAR'
823       include 'COMMON.INTERACT'
824       include 'COMMON.LOCAL'
825       include 'COMMON.NAMES'
826       include 'COMMON.CHAIN'
827       include 'COMMON.FFIELD'
828       include 'COMMON.SBRIDGE'
829       include 'COMMON.HEADER'
830       include 'COMMON.CONTROL'
831 c      include 'COMMON.DBASE'
832 c      include 'COMMON.THREAD'
833       include 'COMMON.TIME1'
834       include 'COMMON.SETUP'
835 C Read bridging residues.
836       read (inp,*) ns,(iss(i),i=1,ns)
837       print *,'ns=',ns
838       if(me.eq.king.or..not.out1file)
839      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
840 C Check whether the specified bridging residues are cystines.
841       do i=1,ns
842         if (itype(iss(i)).ne.1) then
843           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') 
844      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
845      &   ' can form a disulfide bridge?!!!'
846           write (*,'(2a,i3,a)') 
847      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
848      &   ' can form a disulfide bridge?!!!'
849 #ifdef MPI
850          call MPI_Finalize(MPI_COMM_WORLD,ierror)
851          stop
852 #endif
853         endif
854       enddo
855 C Read preformed bridges.
856       if (ns.gt.0) then
857       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
858       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
859       if (nss.gt.0) then
860         nhpb=nss
861 C Check if the residues involved in bridges are in the specified list of
862 C bridging residues.
863         do i=1,nss
864           do j=1,i-1
865             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
866      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
867               write (iout,'(a,i3,a)') 'Disulfide pair',i,
868      &      ' contains residues present in other pairs.'
869               write (*,'(a,i3,a)') 'Disulfide pair',i,
870      &      ' contains residues present in other pairs.'
871 #ifdef MPI
872               call MPI_Finalize(MPI_COMM_WORLD,ierror)
873               stop 
874 #endif
875             endif
876           enddo
877           do j=1,ns
878             if (ihpb(i).eq.iss(j)) goto 10
879           enddo
880           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
881    10     continue
882           do j=1,ns
883             if (jhpb(i).eq.iss(j)) goto 20
884           enddo
885           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
886    20     continue
887           dhpb(i)=dbr
888           forcon(i)=fbr
889         enddo
890         do i=1,nss
891           ihpb(i)=ihpb(i)+nres
892           jhpb(i)=jhpb(i)+nres
893         enddo
894       endif
895       endif
896       return
897       end
898 c----------------------------------------------------------------------------
899       subroutine read_x(kanal,*)
900       implicit real*8 (a-h,o-z)
901       include 'DIMENSIONS'
902       include 'COMMON.GEO'
903       include 'COMMON.VAR'
904       include 'COMMON.CHAIN'
905       include 'COMMON.IOUNITS'
906       include 'COMMON.CONTROL'
907       include 'COMMON.LOCAL'
908       include 'COMMON.INTERACT'
909 c Read coordinates from input
910 c
911       read(kanal,'(8f10.5)',end=10,err=10)
912      &  ((c(l,k),l=1,3),k=1,nres),
913      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
914       do j=1,3
915         c(j,nres+1)=c(j,1)
916         c(j,2*nres)=c(j,nres)
917       enddo
918       call int_from_cart1(.false.)
919       do i=1,nres-1
920         do j=1,3
921           dc(j,i)=c(j,i+1)-c(j,i)
922           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
923         enddo
924       enddo
925       do i=nnt,nct
926         if (itype(i).ne.10) then
927           do j=1,3
928             dc(j,i+nres)=c(j,i+nres)-c(j,i)
929             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
930           enddo
931         endif
932       enddo
933
934       return
935    10 return1
936       end
937 c------------------------------------------------------------------------------
938       subroutine setup_var
939       implicit real*8 (a-h,o-z)
940       include 'DIMENSIONS'
941       include 'COMMON.IOUNITS'
942       include 'COMMON.GEO'
943       include 'COMMON.VAR'
944       include 'COMMON.INTERACT'
945       include 'COMMON.LOCAL'
946       include 'COMMON.NAMES'
947       include 'COMMON.CHAIN'
948       include 'COMMON.FFIELD'
949       include 'COMMON.SBRIDGE'
950       include 'COMMON.HEADER'
951       include 'COMMON.CONTROL'
952 c      include 'COMMON.DBASE'
953 c      include 'COMMON.THREAD'
954       include 'COMMON.TIME1'
955 C Set up variable list.
956       ntheta=nres-2
957       nphi=nres-3
958       nvar=ntheta+nphi
959       nside=0
960       do i=2,nres-1
961         if (itype(i).ne.10) then
962           nside=nside+1
963           ialph(i,1)=nvar+nside
964           ialph(nside,2)=i
965         endif
966       enddo
967       if (indphi.gt.0) then
968         nvar=nphi
969       else if (indback.gt.0) then
970         nvar=nphi+ntheta
971       else
972         nvar=nvar+2*nside
973       endif
974 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
975       return
976       end
977 c----------------------------------------------------------------------------
978       subroutine gen_dist_constr
979 C Generate CA distance constraints.
980       implicit real*8 (a-h,o-z)
981       include 'DIMENSIONS'
982       include 'COMMON.IOUNITS'
983       include 'COMMON.GEO'
984       include 'COMMON.VAR'
985       include 'COMMON.INTERACT'
986       include 'COMMON.LOCAL'
987       include 'COMMON.NAMES'
988       include 'COMMON.CHAIN'
989       include 'COMMON.FFIELD'
990       include 'COMMON.SBRIDGE'
991       include 'COMMON.HEADER'
992       include 'COMMON.CONTROL'
993 c      include 'COMMON.DBASE'
994 c      include 'COMMON.THREAD'
995       include 'COMMON.TIME1'
996       dimension itype_pdb(maxres)
997       common /pizda/ itype_pdb
998       character*2 iden
999 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1000 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1001 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1002 cd     & ' nsup',nsup
1003       do i=nstart_sup,nstart_sup+nsup-1
1004 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1005 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1006         do j=i+2,nstart_sup+nsup-1
1007           nhpb=nhpb+1
1008           ihpb(nhpb)=i+nstart_seq-nstart_sup
1009           jhpb(nhpb)=j+nstart_seq-nstart_sup
1010           forcon(nhpb)=weidis
1011           dhpb(nhpb)=dist(i,j)
1012         enddo
1013       enddo 
1014 cd      write (iout,'(a)') 'Distance constraints:' 
1015 cd      do i=nss+1,nhpb
1016 cd        ii=ihpb(i)
1017 cd        jj=jhpb(i)
1018 cd        iden='CA'
1019 cd        if (ii.gt.nres) then
1020 cd          iden='SC'
1021 cd          ii=ii-nres
1022 cd          jj=jj-nres
1023 cd        endif
1024 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') 
1025 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1026 cd     &  dhpb(i),forcon(i)
1027 cd      enddo
1028       return
1029       end
1030
1031       subroutine read_minim
1032       implicit real*8 (a-h,o-z)
1033       include 'DIMENSIONS'
1034       include 'COMMON.MINIM'
1035       include 'COMMON.IOUNITS'
1036       character*80 ucase
1037       character*320 minimcard
1038       call card_concat(minimcard)
1039       call readi(minimcard,'MAXMIN',maxmin,2000)
1040       call readi(minimcard,'MAXFUN',maxfun,5000)
1041       call readi(minimcard,'MINMIN',minmin,maxmin)
1042       call readi(minimcard,'MINFUN',minfun,maxmin)
1043       call reada(minimcard,'TOLF',tolf,1.0D-2)
1044       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1045       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1046       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1047       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1048       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1049      &         'Options in energy minimization:'
1050       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1051      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1052      & 'MinMin:',MinMin,' MinFun:',MinFun,
1053      & ' TolF:',TolF,' RTolF:',RTolF
1054       return
1055       end
1056 c----------------------------------------------------------------------------
1057       subroutine read_angles(kanal,*)
1058       implicit real*8 (a-h,o-z)
1059       include 'DIMENSIONS'
1060       include 'COMMON.GEO'
1061       include 'COMMON.VAR'
1062       include 'COMMON.CHAIN'
1063       include 'COMMON.IOUNITS'
1064       include 'COMMON.CONTROL'
1065 c Read angles from input 
1066 c
1067        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1068        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1069        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1070        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1071
1072        do i=1,nres
1073 c 9/7/01 avoid 180 deg valence angle
1074         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1075 c
1076         theta(i)=deg2rad*theta(i)
1077         phi(i)=deg2rad*phi(i)
1078         alph(i)=deg2rad*alph(i)
1079         omeg(i)=deg2rad*omeg(i)
1080        enddo
1081       return
1082    10 return1
1083       end
1084 c----------------------------------------------------------------------------
1085       subroutine reada(rekord,lancuch,wartosc,default)
1086       implicit none
1087       character*(*) rekord,lancuch
1088       double precision wartosc,default
1089       integer ilen,iread
1090       external ilen
1091       iread=index(rekord,lancuch)
1092       if (iread.eq.0) then
1093         wartosc=default 
1094         return
1095       endif   
1096       iread=iread+ilen(lancuch)+1
1097       read (rekord(iread:),*,err=10,end=10) wartosc
1098       return
1099   10  wartosc=default
1100       return
1101       end
1102 c----------------------------------------------------------------------------
1103       subroutine readi(rekord,lancuch,wartosc,default)
1104       implicit none
1105       character*(*) rekord,lancuch
1106       integer wartosc,default
1107       integer ilen,iread
1108       external ilen
1109       iread=index(rekord,lancuch)
1110       if (iread.eq.0) then
1111         wartosc=default 
1112         return
1113       endif   
1114       iread=iread+ilen(lancuch)+1
1115       read (rekord(iread:),*,err=10,end=10) wartosc
1116       return
1117   10  wartosc=default
1118       return
1119       end
1120 c----------------------------------------------------------------------------
1121       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1122       implicit none
1123       integer dim,i
1124       integer tablica(dim),default
1125       character*(*) rekord,lancuch
1126       character*80 aux
1127       integer ilen,iread
1128       external ilen
1129       do i=1,dim
1130         tablica(i)=default
1131       enddo
1132       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1133       if (iread.eq.0) return
1134       iread=iread+ilen(lancuch)+1
1135       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1136    10 return
1137       end
1138 c----------------------------------------------------------------------------
1139       subroutine multreada(rekord,lancuch,tablica,dim,default)
1140       implicit none
1141       integer dim,i
1142       double precision tablica(dim),default
1143       character*(*) rekord,lancuch
1144       character*80 aux
1145       integer ilen,iread
1146       external ilen
1147       do i=1,dim
1148         tablica(i)=default
1149       enddo
1150       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1151       if (iread.eq.0) return
1152       iread=iread+ilen(lancuch)+1
1153       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1154    10 return
1155       end
1156 c----------------------------------------------------------------------------
1157       subroutine openunits
1158       implicit real*8 (a-h,o-z)
1159       include 'DIMENSIONS'    
1160 #ifdef MPI
1161       include 'mpif.h'
1162       character*16 form,nodename
1163       integer nodelen
1164 #endif
1165       include 'COMMON.SETUP'
1166       include 'COMMON.IOUNITS'
1167 c      include 'COMMON.MD'
1168       include 'COMMON.CONTROL'
1169       integer lenpre,lenpot,ilen,lentmp
1170       external ilen
1171       character*3 out1file_text,ucase
1172       character*3 ll
1173       external ucase
1174 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1175       call getenv_loc("PREFIX",prefix)
1176       pref_orig = prefix
1177       call getenv_loc("POT",pot)
1178       call getenv_loc("DIRTMP",tmpdir)
1179       call getenv_loc("CURDIR",curdir)
1180       call getenv_loc("OUT1FILE",out1file_text)
1181 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1182       out1file_text=ucase(out1file_text)
1183       if (out1file_text(1:1).eq."Y") then
1184         out1file=.true.
1185       else 
1186         out1file=fg_rank.gt.0
1187       endif
1188       lenpre=ilen(prefix)
1189       lenpot=ilen(pot)
1190       lentmp=ilen(tmpdir)
1191       if (lentmp.gt.0) then
1192           write (*,'(80(1h!))')
1193           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1194           write (*,'(80(1h!))')
1195           write (*,*)"All output files will be on node /tmp directory." 
1196 #ifdef MPI
1197         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1198         if (me.eq.king) then
1199           write (*,*) "The master node is ",nodename
1200         else if (fg_rank.eq.0) then
1201           write (*,*) "I am the CG slave node ",nodename
1202         else 
1203           write (*,*) "I am the FG slave node ",nodename
1204         endif
1205 #endif
1206         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1207         lenpre = lentmp+lenpre+1
1208       endif
1209       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1210 C Get the names and open the input files
1211 #if defined(WINIFL) || defined(WINPGI)
1212       open(1,file=pref_orig(:ilen(pref_orig))//
1213      &  '.inp',status='old',readonly,shared)
1214        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1215 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1216 C Get parameter filenames and open the parameter files.
1217       call getenv_loc('BONDPAR',bondname)
1218       open (ibond,file=bondname,status='old',readonly,shared)
1219       call getenv_loc('THETPAR',thetname)
1220       open (ithep,file=thetname,status='old',readonly,shared)
1221 #ifndef CRYST_THETA
1222       call getenv_loc('THETPARPDB',thetname_pdb)
1223       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1224 #endif
1225       call getenv_loc('ROTPAR',rotname)
1226       open (irotam,file=rotname,status='old',readonly,shared)
1227 #ifndef CRYST_SC
1228       call getenv_loc('ROTPARPDB',rotname_pdb)
1229       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1230 #endif
1231       call getenv_loc('TORPAR',torname)
1232       open (itorp,file=torname,status='old',readonly,shared)
1233       call getenv_loc('TORDPAR',tordname)
1234       open (itordp,file=tordname,status='old',readonly,shared)
1235       call getenv_loc('FOURIER',fouriername)
1236       open (ifourier,file=fouriername,status='old',readonly,shared)
1237       call getenv_loc('ELEPAR',elename)
1238       open (ielep,file=elename,status='old',readonly,shared)
1239       call getenv_loc('SIDEPAR',sidename)
1240       open (isidep,file=sidename,status='old',readonly,shared)
1241 #elif (defined CRAY) || (defined AIX)
1242       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1243      &  action='read')
1244 c      print *,"Processor",myrank," opened file 1" 
1245       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1246 c      print *,"Processor",myrank," opened file 9" 
1247 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1248 C Get parameter filenames and open the parameter files.
1249       call getenv_loc('BONDPAR',bondname)
1250       open (ibond,file=bondname,status='old',action='read')
1251 c      print *,"Processor",myrank," opened file IBOND" 
1252       call getenv_loc('THETPAR',thetname)
1253       open (ithep,file=thetname,status='old',action='read')
1254 c      print *,"Processor",myrank," opened file ITHEP" 
1255 #ifndef CRYST_THETA
1256       call getenv_loc('THETPARPDB',thetname_pdb)
1257       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1258 #endif
1259       call getenv_loc('ROTPAR',rotname)
1260       open (irotam,file=rotname,status='old',action='read')
1261 c      print *,"Processor",myrank," opened file IROTAM" 
1262 #ifndef CRYST_SC
1263       call getenv_loc('ROTPARPDB',rotname_pdb)
1264       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1265 #endif
1266       call getenv_loc('TORPAR',torname)
1267       open (itorp,file=torname,status='old',action='read')
1268 c      print *,"Processor",myrank," opened file ITORP" 
1269       call getenv_loc('TORDPAR',tordname)
1270       open (itordp,file=tordname,status='old',action='read')
1271 c      print *,"Processor",myrank," opened file ITORDP" 
1272       call getenv_loc('SCCORPAR',sccorname)
1273       open (isccor,file=sccorname,status='old',action='read')
1274 c      print *,"Processor",myrank," opened file ISCCOR" 
1275       call getenv_loc('FOURIER',fouriername)
1276       open (ifourier,file=fouriername,status='old',action='read')
1277 c      print *,"Processor",myrank," opened file IFOURIER" 
1278       call getenv_loc('ELEPAR',elename)
1279       open (ielep,file=elename,status='old',action='read')
1280 c      print *,"Processor",myrank," opened file IELEP" 
1281       call getenv_loc('SIDEPAR',sidename)
1282       open (isidep,file=sidename,status='old',action='read')
1283 c      print *,"Processor",myrank," opened file ISIDEP" 
1284 c      print *,"Processor",myrank," opened parameter files" 
1285 #elif (defined G77)
1286       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1287       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1288 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1289 C Get parameter filenames and open the parameter files.
1290       call getenv_loc('BONDPAR',bondname)
1291       open (ibond,file=bondname,status='old')
1292       call getenv_loc('THETPAR',thetname)
1293       open (ithep,file=thetname,status='old')
1294 #ifndef CRYST_THETA
1295       call getenv_loc('THETPARPDB',thetname_pdb)
1296       open (ithep_pdb,file=thetname_pdb,status='old')
1297 #endif
1298       call getenv_loc('ROTPAR',rotname)
1299       open (irotam,file=rotname,status='old')
1300 #ifndef CRYST_SC
1301       call getenv_loc('ROTPARPDB',rotname_pdb)
1302       open (irotam_pdb,file=rotname_pdb,status='old')
1303 #endif
1304       call getenv_loc('TORPAR',torname)
1305       open (itorp,file=torname,status='old')
1306       call getenv_loc('TORDPAR',tordname)
1307       open (itordp,file=tordname,status='old')
1308       call getenv_loc('SCCORPAR',sccorname)
1309       open (isccor,file=sccorname,status='old')
1310       call getenv_loc('FOURIER',fouriername)
1311       open (ifourier,file=fouriername,status='old')
1312       call getenv_loc('ELEPAR',elename)
1313       open (ielep,file=elename,status='old')
1314       call getenv_loc('SIDEPAR',sidename)
1315       open (isidep,file=sidename,status='old')
1316 #else
1317       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1318      &  readonly)
1319        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1320 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1321 C Get parameter filenames and open the parameter files.
1322       call getenv_loc('BONDPAR',bondname)
1323       open (ibond,file=bondname,status='old',readonly)
1324       call getenv_loc('THETPAR',thetname)
1325       open (ithep,file=thetname,status='old',readonly)
1326 #ifndef CRYST_THETA
1327       call getenv_loc('THETPARPDB',thetname_pdb)
1328       print *,"thetname_pdb ",thetname_pdb
1329       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1330       print *,ithep_pdb," opened"
1331 #endif
1332       call getenv_loc('ROTPAR',rotname)
1333       open (irotam,file=rotname,status='old',readonly)
1334 #ifndef CRYST_SC
1335       call getenv_loc('ROTPARPDB',rotname_pdb)
1336       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1337 #endif
1338       call getenv_loc('TORPAR',torname)
1339       open (itorp,file=torname,status='old',readonly)
1340       call getenv_loc('TORDPAR',tordname)
1341       open (itordp,file=tordname,status='old',readonly)
1342       call getenv_loc('SCCORPAR',sccorname)
1343       open (isccor,file=sccorname,status='old',readonly)
1344       call getenv_loc('FOURIER',fouriername)
1345       open (ifourier,file=fouriername,status='old',readonly)
1346       call getenv_loc('ELEPAR',elename)
1347       open (ielep,file=elename,status='old',readonly)
1348       call getenv_loc('SIDEPAR',sidename)
1349       open (isidep,file=sidename,status='old',readonly)
1350 #endif
1351 #ifndef OLDSCP
1352 C
1353 C 8/9/01 In the newest version SCp interaction constants are read from a file
1354 C Use -DOLDSCP to use hard-coded constants instead.
1355 C
1356       call getenv_loc('SCPPAR',scpname)
1357 #if defined(WINIFL) || defined(WINPGI)
1358       open (iscpp,file=scpname,status='old',readonly,shared)
1359 #elif (defined CRAY)  || (defined AIX)
1360       open (iscpp,file=scpname,status='old',action='read')
1361 #elif (defined G77)
1362       open (iscpp,file=scpname,status='old')
1363 #else
1364       open (iscpp,file=scpname,status='old',readonly)
1365 #endif
1366 #endif
1367       call getenv_loc('PATTERN',patname)
1368 #if defined(WINIFL) || defined(WINPGI)
1369       open (icbase,file=patname,status='old',readonly,shared)
1370 #elif (defined CRAY)  || (defined AIX)
1371       open (icbase,file=patname,status='old',action='read')
1372 #elif (defined G77)
1373       open (icbase,file=patname,status='old')
1374 #else
1375       open (icbase,file=patname,status='old',readonly)
1376 #endif
1377 #ifdef MPI
1378 C Open output file only for CG processes
1379 c      print *,"Processor",myrank," fg_rank",fg_rank
1380       if (fg_rank.eq.0) then
1381
1382       if (nodes.eq.1) then
1383         npos=3
1384       else
1385         npos = dlog10(dfloat(nodes-1))+1
1386       endif
1387       if (npos.lt.3) npos=3
1388       write (liczba,'(i1)') npos
1389       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1390      &  //')'
1391       write (liczba,form) me
1392       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1393      &  liczba(:ilen(liczba))
1394       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1395      &  //'.int'
1396       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1397      &  //'.pdb'
1398       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1399      &  liczba(:ilen(liczba))//'.mol2'
1400       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1401      &  liczba(:ilen(liczba))//'.stat'
1402       if (lentmp.gt.0)
1403      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1404      &      //liczba(:ilen(liczba))//'.stat')
1405       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1406      &  //'.rst'
1407 c      if(usampl) then
1408 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1409 c     & liczba(:ilen(liczba))//'.const'
1410 c      endif 
1411
1412       endif
1413 #else
1414       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1415       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1416       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1417       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1418       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1419       if (lentmp.gt.0)
1420      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1421      &    '.stat')
1422       rest2name=prefix(:ilen(prefix))//'.rst'
1423 c      if(usampl) then 
1424 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1425 c      endif 
1426 #endif
1427 #if defined(AIX) || defined(PGI)
1428       if (me.eq.king .or. .not. out1file) 
1429      &   open(iout,file=outname,status='unknown')
1430 #ifdef DEBUG
1431       if (fg_rank.gt.0) then
1432         write (liczba,'(i3.3)') myrank/nfgtasks
1433         write (ll,'(bz,i3.3)') fg_rank
1434         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1435      &   status='unknown')
1436       endif
1437 #endif
1438       if(me.eq.king) then
1439        open(igeom,file=intname,status='unknown',position='append')
1440        open(ipdb,file=pdbname,status='unknown')
1441        open(imol2,file=mol2name,status='unknown')
1442        open(istat,file=statname,status='unknown',position='append')
1443       else
1444 c1out       open(iout,file=outname,status='unknown')
1445       endif
1446 #else
1447       if (me.eq.king .or. .not.out1file)
1448      &    open(iout,file=outname,status='unknown')
1449 #ifdef DEBUG
1450       if (fg_rank.gt.0) then
1451         write (liczba,'(i3.3)') myrank/nfgtasks
1452         write (ll,'(bz,i3.3)') fg_rank
1453         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1454      &   status='unknown')
1455       endif
1456 #endif
1457       if(me.eq.king) then
1458        open(igeom,file=intname,status='unknown',access='append')
1459        open(ipdb,file=pdbname,status='unknown')
1460        open(imol2,file=mol2name,status='unknown')
1461        open(istat,file=statname,status='unknown',access='append')
1462       else
1463 c1out       open(iout,file=outname,status='unknown')
1464       endif
1465 #endif
1466       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1467       csa_seed=prefix(:lenpre)//'.CSA.seed'
1468       csa_history=prefix(:lenpre)//'.CSA.history'
1469       csa_bank=prefix(:lenpre)//'.CSA.bank'
1470       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1471       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1472       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1473 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1474       csa_int=prefix(:lenpre)//'.int'
1475       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1476       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1477       csa_in=prefix(:lenpre)//'.CSA.in'
1478 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1479 C Write file names
1480       if (me.eq.king)then
1481       write (iout,'(80(1h-))')
1482       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1483       write (iout,'(80(1h-))')
1484       write (iout,*) "Input file                      : ",
1485      &  pref_orig(:ilen(pref_orig))//'.inp'
1486       write (iout,*) "Output file                     : ",
1487      &  outname(:ilen(outname))
1488       write (iout,*)
1489       write (iout,*) "Sidechain potential file        : ",
1490      &  sidename(:ilen(sidename))
1491 #ifndef OLDSCP
1492       write (iout,*) "SCp potential file              : ",
1493      &  scpname(:ilen(scpname))
1494 #endif
1495       write (iout,*) "Electrostatic potential file    : ",
1496      &  elename(:ilen(elename))
1497       write (iout,*) "Cumulant coefficient file       : ",
1498      &  fouriername(:ilen(fouriername))
1499       write (iout,*) "Torsional parameter file        : ",
1500      &  torname(:ilen(torname))
1501       write (iout,*) "Double torsional parameter file : ",
1502      &  tordname(:ilen(tordname))
1503       write (iout,*) "SCCOR parameter file : ",
1504      &  sccorname(:ilen(sccorname))
1505       write (iout,*) "Bond & inertia constant file    : ",
1506      &  bondname(:ilen(bondname))
1507       write (iout,*) "Bending parameter file          : ",
1508      &  thetname(:ilen(thetname))
1509       write (iout,*) "Rotamer parameter file          : ",
1510      &  rotname(:ilen(rotname))
1511       write (iout,*) "Threading database              : ",
1512      &  patname(:ilen(patname))
1513       if (lentmp.ne.0) 
1514      &write (iout,*)" DIRTMP                          : ",
1515      &  tmpdir(:lentmp)
1516       write (iout,'(80(1h-))')
1517       endif
1518       return
1519       end
1520 c----------------------------------------------------------------------------
1521       subroutine card_concat(card)
1522       implicit real*8 (a-h,o-z)
1523       include 'DIMENSIONS'
1524       include 'COMMON.IOUNITS'
1525       character*(*) card
1526       character*80 karta,ucase
1527       external ilen
1528       read (inp,'(a)') karta
1529       karta=ucase(karta)
1530       card=' '
1531       do while (karta(80:80).eq.'&')
1532         card=card(:ilen(card)+1)//karta(:79)
1533         read (inp,'(a)') karta
1534         karta=ucase(karta)
1535       enddo
1536       card=card(:ilen(card)+1)//karta
1537       return
1538       end
1539
1540       subroutine read_dist_constr
1541       implicit real*8 (a-h,o-z)
1542       include 'DIMENSIONS'
1543 #ifdef MPI
1544       include 'mpif.h'
1545 #endif
1546       include 'COMMON.SETUP'
1547       include 'COMMON.CONTROL'
1548       include 'COMMON.CHAIN'
1549       include 'COMMON.IOUNITS'
1550       include 'COMMON.SBRIDGE'
1551       integer ifrag_(2,100),ipair_(2,100)
1552       double precision wfrag_(100),wpair_(100)
1553       character*500 controlcard
1554 c      write (iout,*) "Calling read_dist_constr"
1555 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1556 c      call flush(iout)
1557       call card_concat(controlcard)
1558       call readi(controlcard,"NFRAG",nfrag_,0)
1559       call readi(controlcard,"NPAIR",npair_,0)
1560       call readi(controlcard,"NDIST",ndist_,0)
1561       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1562       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1563       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1564       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1565       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1566 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1567 c      write (iout,*) "IFRAG"
1568 c      do i=1,nfrag_
1569 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1570 c      enddo
1571 c      write (iout,*) "IPAIR"
1572 c      do i=1,npair_
1573 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1574 c      enddo
1575       call flush(iout)
1576       do i=1,nfrag_
1577         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1578         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1579      &    ifrag_(2,i)=nstart_sup+nsup-1
1580 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1581         call flush(iout)
1582         if (wfrag_(i).gt.0.0d0) then
1583         do j=ifrag_(1,i),ifrag_(2,i)-1
1584           do k=j+1,ifrag_(2,i)
1585             write (iout,*) "j",j," k",k
1586             ddjk=dist(j,k)
1587             if (constr_dist.eq.1) then
1588             nhpb=nhpb+1
1589             ihpb(nhpb)=j
1590             jhpb(nhpb)=k
1591               dhpb(nhpb)=ddjk
1592             forcon(nhpb)=wfrag_(i) 
1593             else if (constr_dist.eq.2) then
1594               if (ddjk.le.dist_cut) then
1595                 nhpb=nhpb+1
1596                 ihpb(nhpb)=j
1597                 jhpb(nhpb)=k
1598                 dhpb(nhpb)=ddjk
1599                 forcon(nhpb)=wfrag_(i) 
1600               endif
1601             else
1602               nhpb=nhpb+1
1603               ihpb(nhpb)=j
1604               jhpb(nhpb)=k
1605               dhpb(nhpb)=ddjk
1606               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1607             endif
1608 #ifdef MPI
1609             if (.not.out1file .or. me.eq.king) 
1610      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1611      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1612 #else
1613             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1614      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1615 #endif
1616           enddo
1617         enddo
1618         endif
1619       enddo
1620       do i=1,npair_
1621         if (wpair_(i).gt.0.0d0) then
1622         ii = ipair_(1,i)
1623         jj = ipair_(2,i)
1624         if (ii.gt.jj) then
1625           itemp=ii
1626           ii=jj
1627           jj=itemp
1628         endif
1629         do j=ifrag_(1,ii),ifrag_(2,ii)
1630           do k=ifrag_(1,jj),ifrag_(2,jj)
1631             nhpb=nhpb+1
1632             ihpb(nhpb)=j
1633             jhpb(nhpb)=k
1634             forcon(nhpb)=wpair_(i)
1635             dhpb(nhpb)=dist(j,k)
1636 #ifdef MPI
1637             if (.not.out1file .or. me.eq.king)
1638      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1639      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1640 #else
1641             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1642      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1643 #endif
1644           enddo
1645         enddo
1646         endif
1647       enddo 
1648       do i=1,ndist_
1649         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1650         if (forcon(nhpb+1).gt.0.0d0) then
1651           nhpb=nhpb+1
1652           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1653 #ifdef MPI
1654           if (.not.out1file .or. me.eq.king)
1655      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1656      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1657 #else
1658           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1659      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1660 #endif
1661         endif
1662       enddo
1663       call flush(iout)
1664       return
1665       end
1666 c-------------------------------------------------------------------------------
1667 #ifdef WINIFL
1668       subroutine flush(iu)
1669       return
1670       end
1671 #endif
1672 #ifdef AIX
1673       subroutine flush(iu)
1674       call flush_(iu)
1675       return
1676       end
1677 #endif
1678 c------------------------------------------------------------------------------
1679       subroutine copy_to_tmp(source)
1680       include "DIMENSIONS"
1681       include "COMMON.IOUNITS"
1682       character*(*) source
1683       character* 256 tmpfile
1684       integer ilen
1685       external ilen
1686       logical ex
1687       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1688       inquire(file=tmpfile,exist=ex)
1689       if (ex) then
1690         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1691      &   " to temporary directory..."
1692         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1693         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1694       endif
1695       return
1696       end
1697 c------------------------------------------------------------------------------
1698       subroutine move_from_tmp(source)
1699       include "DIMENSIONS"
1700       include "COMMON.IOUNITS"
1701       character*(*) source
1702       integer ilen
1703       external ilen
1704       write (*,*) "Moving ",source(:ilen(source)),
1705      & " from temporary directory to working directory"
1706       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1707       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1708       return
1709       end
1710 c------------------------------------------------------------------------------
1711       subroutine random_init(seed)
1712 C
1713 C Initialize random number generator
1714 C
1715       implicit real*8 (a-h,o-z)
1716       include 'DIMENSIONS'
1717 #ifdef AMD64
1718       integer*8 iseedi8
1719 #endif
1720 #ifdef MPI
1721       include 'mpif.h'
1722       logical OKRandom, prng_restart
1723       real*8  r1
1724       integer iseed_array(4)
1725 #endif
1726       include 'COMMON.IOUNITS'
1727       include 'COMMON.TIME1'
1728 c      include 'COMMON.THREAD'
1729       include 'COMMON.SBRIDGE'
1730       include 'COMMON.CONTROL'
1731       include 'COMMON.MCM'
1732 c      include 'COMMON.MAP'
1733       include 'COMMON.HEADER'
1734 c      include 'COMMON.CSA'
1735       include 'COMMON.CHAIN'
1736 c      include 'COMMON.MUCA'
1737 c      include 'COMMON.MD'
1738       include 'COMMON.FFIELD'
1739       include 'COMMON.SETUP'
1740       iseed=-dint(dabs(seed))
1741       if (iseed.eq.0) then
1742         write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1743      &    'Random seed undefined. The program will stop.'
1744         write (*,'(/80(1h*)/20x,a/80(1h*))') 
1745      &    'Random seed undefined. The program will stop.'
1746 #ifdef MPI
1747         call mpi_finalize(mpi_comm_world,ierr)
1748 #endif
1749         stop 'Bad random seed.'
1750       endif
1751 #ifdef MPI
1752       if (fg_rank.eq.0) then
1753       seed=seed*(me+1)+1
1754 #ifdef AMD64
1755       iseedi8=dint(seed)
1756       if(me.eq.king .or. .not. out1file)
1757      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1758       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1759       OKRandom = prng_restart(me,iseedi8)
1760 #else
1761       do i=1,4
1762        tmp=65536.0d0**(4-i)
1763        iseed_array(i) = dint(seed/tmp)
1764        seed=seed-iseed_array(i)*tmp
1765       enddo
1766       if(me.eq.king .or. .not. out1file)
1767      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1768      &                 (iseed_array(i),i=1,4)
1769       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1770      &                 (iseed_array(i),i=1,4)
1771       OKRandom = prng_restart(me,iseed_array)
1772 #endif
1773       if (OKRandom) then
1774 c        r1=ran_number(0.0D0,1.0D0)
1775         if(me.eq.king .or. .not. out1file)
1776      &   write (iout,*) 'ran_num',r1
1777         if (r1.lt.0.0d0) OKRandom=.false.
1778       endif
1779       if (.not.OKRandom) then
1780         write (iout,*) 'PRNG IS NOT WORKING!!!'
1781         print *,'PRNG IS NOT WORKING!!!'
1782         if (me.eq.0) then 
1783          call flush(iout)
1784          call mpi_abort(mpi_comm_world,error_msg,ierr)
1785          stop
1786         else
1787          write (iout,*) 'too many processors for parallel prng'
1788          write (*,*) 'too many processors for parallel prng'
1789          call flush(iout)
1790          stop
1791         endif
1792       endif
1793       endif
1794 #else
1795       call vrndst(iseed)
1796 c      write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1797 #endif
1798       return
1799       end