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