68ae6f415ea08d1293184b9a695bd226814cc8a6
[unres.git] / readrtns_CSA.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       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21       if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23       if (modecalc.eq.14) then 
24          call read_MDpar
25          call read_REMDpar
26       endif
27 C Read MUCA control parameters if required
28       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       if (modecalc.eq.8) then
32        inquire (file="fort.40",exist=file_exist)
33        if (.not.file_exist) call csaread
34       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.THREAD'
70       include 'COMMON.SBRIDGE'
71       include 'COMMON.CONTROL'
72       include 'COMMON.MCM'
73       include 'COMMON.MAP'
74       include 'COMMON.HEADER'
75       include 'COMMON.CSA'
76       include 'COMMON.CHAIN'
77       include 'COMMON.MUCA'
78       include 'COMMON.MD'
79       include 'COMMON.FFIELD'
80       include 'COMMON.INTERACT'
81       include 'COMMON.SETUP'
82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
84       character*80 ucase
85       character*320 controlcard
86
87       nglob_csa=0
88       eglob_csa=1d99
89       nmin_csa=0
90       read (INP,'(a)') titel
91       call card_concat(controlcard)
92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94       call reada(controlcard,'SEED',seed,0.0D0)
95       call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97       read_cart=index(controlcard,'READ_CART').gt.0
98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99       call readi(controlcard,'SYM',symetr,1)
100       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107       call reada(controlcard,'DRMS',drms,0.1D0)
108       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc 
110        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 
111        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max 
112        write (iout,'(a,f10.1)')'DRMS    = ',drms 
113        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm 
114        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
115       endif
116       call readi(controlcard,'NZ_START',nz_start,0)
117       call readi(controlcard,'NZ_END',nz_end,0)
118       call readi(controlcard,'IZ_SC',iz_sc,0)
119       timlim=60.0D0*timlim
120       safety = 60.0d0*safety
121       timem=timlim
122       modecalc=0
123       call reada(controlcard,"T_BATH",t_bath,300.0d0)
124       minim=(index(controlcard,'MINIMIZE').gt.0)
125       dccart=(index(controlcard,'CART').gt.0)
126       overlapsc=(index(controlcard,'OVERLAP').gt.0)
127       overlapsc=.not.overlapsc
128       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129       searchsc=.not.searchsc
130       sideadd=(index(controlcard,'SIDEADD').gt.0)
131       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132       outpdb=(index(controlcard,'PDBOUT').gt.0)
133       outmol2=(index(controlcard,'MOL2OUT').gt.0)
134       pdbref=(index(controlcard,'PDBREF').gt.0)
135       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136       indpdb=index(controlcard,'PDBSTART')
137       extconf=(index(controlcard,'EXTCONF').gt.0)
138       call readi(controlcard,'IPRINT',iprint,0)
139       call readi(controlcard,'MAXGEN',maxgen,10000)
140       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141       call readi(controlcard,"KDIAG",kdiag,0)
142       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144      & write (iout,*) "RESCALE_MODE",rescale_mode
145       split_ene=index(controlcard,'SPLIT_ENE').gt.0
146       if (index(controlcard,'REGULAR').gt.0.0D0) then
147         call reada(controlcard,'WEIDIS',weidis,0.1D0)
148         modecalc=1
149         refstr=.true.
150       endif
151       if (index(controlcard,'CHECKGRAD').gt.0) then
152         modecalc=5
153         if (index(controlcard,'CART').gt.0) then
154           icheckgrad=1
155         elseif (index(controlcard,'CARINT').gt.0) then
156           icheckgrad=2
157         else
158           icheckgrad=3
159         endif
160       elseif (index(controlcard,'THREAD').gt.0) then
161         modecalc=2
162         call readi(controlcard,'THREAD',nthread,0)
163         if (nthread.gt.0) then
164           call reada(controlcard,'WEIDIS',weidis,0.1D0)
165         else
166           if (fg_rank.eq.0)
167      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
168           stop 'Error termination in Read_Control.'
169         endif
170       else if (index(controlcard,'MCMA').gt.0) then
171         modecalc=3
172       else if (index(controlcard,'MCEE').gt.0) then
173         modecalc=6
174       else if (index(controlcard,'MULTCONF').gt.0) then
175         modecalc=4
176       else if (index(controlcard,'MAP').gt.0) then
177         modecalc=7
178         call readi(controlcard,'MAP',nmap,0)
179       else if (index(controlcard,'CSA').gt.0) then
180         modecalc=8
181 crc      else if (index(controlcard,'ZSCORE').gt.0) then
182 crc   
183 crc  ZSCORE is rm from UNRES, modecalc=9 is available
184 crc
185 crc        modecalc=9
186 cfcm      else if (index(controlcard,'MCMF').gt.0) then
187 cfmc        modecalc=10
188       else if (index(controlcard,'SOFTREG').gt.0) then
189         modecalc=11
190       else if (index(controlcard,'CHECK_BOND').gt.0) then
191         modecalc=-1
192       else if (index(controlcard,'TEST').gt.0) then
193         modecalc=-2
194       else if (index(controlcard,'MD').gt.0) then
195         modecalc=12
196       else if (index(controlcard,'RE ').gt.0) then
197         modecalc=14
198       endif
199
200       lmuca=index(controlcard,'MUCA').gt.0
201       call readi(controlcard,'MUCADYN',mucadyn,0)      
202       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203       if (lmuca .and. (me.eq.king .or. .not.out1file )) 
204      & then
205        write (iout,*) 'MUCADYN=',mucadyn
206        write (iout,*) 'MUCASMOOTH=',muca_smooth
207       endif
208
209       iscode=index(controlcard,'ONE_LETTER')
210       indphi=index(controlcard,'PHI')
211       indback=index(controlcard,'BACK')
212       iranconf=index(controlcard,'RAND_CONF')
213       i2ndstr=index(controlcard,'USE_SEC_PRED')
214       gradout=index(controlcard,'GRADOUT').gt.0
215       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217       if (me.eq.king .or. .not.out1file ) 
218      &  write (iout,*) "DISTCHAINMAX",distchainmax
219       
220       if(me.eq.king.or..not.out1file)
221      & write (iout,'(2a)') diagmeth(kdiag),
222      &  ' routine used to diagonalize matrices.'
223       return
224       end
225 c--------------------------------------------------------------------------
226       subroutine read_REMDpar
227 C
228 C Read REMD settings
229 C
230       implicit real*8 (a-h,o-z)
231       include 'DIMENSIONS'
232       include 'COMMON.IOUNITS'
233       include 'COMMON.TIME1'
234       include 'COMMON.MD'
235 #ifndef LANG0
236       include 'COMMON.LANGEVIN'
237 #else
238       include 'COMMON.LANGEVIN.lang0'
239 #endif
240       include 'COMMON.INTERACT'
241       include 'COMMON.NAMES'
242       include 'COMMON.GEO'
243       include 'COMMON.REMD'
244       include 'COMMON.CONTROL'
245       include 'COMMON.SETUP'
246       character*80 ucase
247       character*320 controlcard
248       character*3200 controlcard1
249       integer iremd_m_total
250
251       if(me.eq.king.or..not.out1file)
252      & write (iout,*) "REMD setup"
253
254       call card_concat(controlcard)
255       call readi(controlcard,"NREP",nrep,3)
256       call readi(controlcard,"NSTEX",nstex,1000)
257       call reada(controlcard,"RETMIN",retmin,10.0d0)
258       call reada(controlcard,"RETMAX",retmax,1000.0d0)
259       mremdsync=(index(controlcard,'SYNC').gt.0)
260       call readi(controlcard,"NSYN",i_sync_step,100)
261       restart1file=(index(controlcard,'REST1FILE').gt.0)
262       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264       if(max_cache_traj_use.gt.max_cache_traj)
265      &           max_cache_traj_use=max_cache_traj
266       if(me.eq.king.or..not.out1file) then
267 cd       if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd        write (iout,*) "NTWX value is ignored"
270 cd        write (iout,*) "  trajectory is stored to one file by master"
271 cd        write (iout,*) "  before exchange at NSTEX intervals"
272 cd       endif
273        write (iout,*) "NREP= ",nrep
274        write (iout,*) "NSTEX= ",nstex
275        write (iout,*) "SYNC= ",mremdsync 
276        write (iout,*) "NSYN= ",i_sync_step
277        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
278       endif
279       remd_tlist=.false.
280       if (index(controlcard,'TLIST').gt.0) then
281          remd_tlist=.true.
282          call card_concat(controlcard1)
283          read(controlcard1,*) (remd_t(i),i=1,nrep) 
284          if(me.eq.king.or..not.out1file)
285      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep) 
286       endif
287       remd_mlist=.false.
288       if (index(controlcard,'MLIST').gt.0) then
289          remd_mlist=.true.
290          call card_concat(controlcard1)
291          read(controlcard1,*) (remd_m(i),i=1,nrep)  
292          if(me.eq.king.or..not.out1file) then
293           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
294           iremd_m_total=0
295           do i=1,nrep
296            iremd_m_total=iremd_m_total+remd_m(i)
297           enddo
298           write (iout,*) 'Total number of replicas ',iremd_m_total
299          endif
300       endif
301       if(me.eq.king.or..not.out1file) 
302      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
303       return
304       end
305 c--------------------------------------------------------------------------
306       subroutine read_MDpar
307 C
308 C Read MD settings
309 C
310       implicit real*8 (a-h,o-z)
311       include 'DIMENSIONS'
312       include 'COMMON.IOUNITS'
313       include 'COMMON.TIME1'
314       include 'COMMON.MD'
315 #ifndef LANG0
316       include 'COMMON.LANGEVIN'
317 #else
318       include 'COMMON.LANGEVIN.lang0'
319 #endif
320       include 'COMMON.INTERACT'
321       include 'COMMON.NAMES'
322       include 'COMMON.GEO'
323       include 'COMMON.SETUP'
324       include 'COMMON.CONTROL'
325       include 'COMMON.SPLITELE'
326       character*80 ucase
327       character*320 controlcard
328
329       call card_concat(controlcard)
330       call readi(controlcard,"NSTEP",n_timestep,1000000)
331       call readi(controlcard,"NTWE",ntwe,100)
332       call readi(controlcard,"NTWX",ntwx,1000)
333       call reada(controlcard,"DT",d_time,1.0d-1)
334       call reada(controlcard,"DVMAX",dvmax,2.0d1)
335       call reada(controlcard,"DAMAX",damax,1.0d1)
336       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337       call readi(controlcard,"LANG",lang,0)
338       RESPA = index(controlcard,"RESPA") .gt. 0
339       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340       ntime_split0=ntime_split
341       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342       ntime_split0=ntime_split
343       call reada(controlcard,"R_CUT",r_cut,2.0d0)
344       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345       rest = index(controlcard,"REST").gt.0
346       tbf = index(controlcard,"TBF").gt.0
347       usampl = index(controlcard,"USAMPL").gt.0
348       mdpdb = index(controlcard,"MDPDB").gt.0
349       call reada(controlcard,"T_BATH",t_bath,300.0d0)
350       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) 
351       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
352       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
353       if (count_reset_moment.eq.0) count_reset_moment=1000000000
354       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
355       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
356       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
357       if (count_reset_vel.eq.0) count_reset_vel=1000000000
358       large = index(controlcard,"LARGE").gt.0
359       print_compon = index(controlcard,"PRINT_COMPON").gt.0
360       rattle = index(controlcard,"RATTLE").gt.0
361 c  if performing umbrella sampling, fragments constrained are read from the fragment file 
362       nset=0
363       if(usampl) then
364         call read_fragments
365       endif
366       
367       if(me.eq.king.or..not.out1file) then
368        write (iout,*)
369        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
370        write (iout,*)
371        write (iout,'(a)') "The units are:"
372        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
373        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
374      &  " acceleration: angstrom/(48.9 fs)**2"
375        write (iout,'(a)') "energy: kcal/mol, temperature: K"
376        write (iout,*)
377        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
378        write (iout,'(a60,f10.5,a)') 
379      &  "Initial time step of numerical integration:",d_time,
380      &  " natural units"
381        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
382        if (RESPA) then
383         write (iout,'(2a,i4,a)') 
384      &    "A-MTS algorithm used; initial time step for fast-varying",
385      &    " short-range forces split into",ntime_split," steps."
386         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
387      &   r_cut," lambda",rlamb
388        endif
389        write (iout,'(2a,f10.5)') 
390      &  "Maximum acceleration threshold to reduce the time step",
391      &  "/increase split number:",damax
392        write (iout,'(2a,f10.5)') 
393      &  "Maximum predicted energy drift to reduce the timestep",
394      &  "/increase split number:",edriftmax
395        write (iout,'(a60,f10.5)') 
396      & "Maximum velocity threshold to reduce velocities:",dvmax
397        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
398        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
399        if (rattle) write (iout,'(a60)') 
400      &  "Rattle algorithm used to constrain the virtual bonds"
401       endif
402       reset_fricmat=1000
403       if (lang.gt.0) then
404         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
405         call reada(controlcard,"RWAT",rwat,1.4d0)
406         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
407         surfarea=index(controlcard,"SURFAREA").gt.0
408         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
409         if(me.eq.king.or..not.out1file)then
410          write (iout,'(/a,$)') "Langevin dynamics calculation"
411          if (lang.eq.1) then
412           write (iout,'(a/)') 
413      &      " with direct integration of Langevin equations"  
414          else if (lang.eq.2) then
415           write (iout,'(a/)') " with TINKER stochasic MD integrator"
416          else if (lang.eq.3) then
417           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
418          else if (lang.eq.4) then
419           write (iout,'(a/)') " in overdamped mode"
420          else
421           write (iout,'(//a,i5)') 
422      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
423           stop
424          endif
425          write (iout,'(a60,f10.5)') "Temperature:",t_bath
426          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
427          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
428          write (iout,'(a60,f10.5)') 
429      &   "Scaling factor of the friction forces:",scal_fric
430          if (surfarea) write (iout,'(2a,i10,a)') 
431      &     "Friction coefficients will be scaled by solvent-accessible",
432      &     " surface area every",reset_fricmat," steps."
433         endif
434 c Calculate friction coefficients and bounds of stochastic forces
435         eta=6*pi*cPoise*etawat
436         if(me.eq.king.or..not.out1file)
437      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
438      &   ,eta
439         gamp=scal_fric*(pstok+rwat)*eta
440         stdfp=dsqrt(2*Rb*t_bath/d_time)
441         do i=1,ntyp
442           gamsc(i)=scal_fric*(restok(i)+rwat)*eta  
443           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
444         enddo 
445         if(me.eq.king.or..not.out1file)then
446          write (iout,'(/2a/)') 
447      &   "Radii of site types and friction coefficients and std's of",
448      &   " stochastic forces of fully exposed sites"
449          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
450          do i=1,ntyp
451           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
452      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
453          enddo
454         endif
455       else if (tbf) then
456         if(me.eq.king.or..not.out1file)then
457          write (iout,'(a)') "Berendsen bath calculation"
458          write (iout,'(a60,f10.5)') "Temperature:",t_bath
459          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
460          if (reset_moment) 
461      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
462      &   count_reset_moment," steps"
463          if (reset_vel) 
464      &    write (iout,'(a,i10,a)') 
465      &    "Velocities will be reset at random every",count_reset_vel,
466      &   " steps"
467         endif
468       else
469         if(me.eq.king.or..not.out1file)
470      &   write (iout,'(a31)') "Microcanonical mode calculation"
471       endif
472       if(me.eq.king.or..not.out1file)then
473        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
474        if (usampl) then
475           write(iout,*) "MD running with constraints."
476           write(iout,*) "Equilibration time ", eq_time, " mtus." 
477           write(iout,*) "Constraining ", nfrag," fragments."
478           write(iout,*) "Length of each fragment, weight and q0:"
479           do iset=1,nset
480            write (iout,*) "Set of restraints #",iset
481            do i=1,nfrag
482               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
483      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
484            enddo
485            write(iout,*) "constraints between ", npair, "fragments."
486            write(iout,*) "constraint pairs, weights and q0:"
487            do i=1,npair
488             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
489      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
490            enddo
491            write(iout,*) "angle constraints within ", nfrag_back, 
492      &      "backbone fragments."
493            write(iout,*) "fragment, weights:"
494            do i=1,nfrag_back
495             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
496      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
497      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
498            enddo
499           enddo
500         iset=mod(kolor,nset)+1
501        endif
502       endif
503       if(me.eq.king.or..not.out1file)
504      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
505       return
506       end
507 c------------------------------------------------------------------------------
508       subroutine molread
509 C
510 C Read molecular data.
511 C
512       implicit real*8 (a-h,o-z)
513       include 'DIMENSIONS'
514 #ifdef MPI
515       include 'mpif.h'
516       integer error_msg
517 #endif
518       include 'COMMON.IOUNITS'
519       include 'COMMON.GEO'
520       include 'COMMON.VAR'
521       include 'COMMON.INTERACT'
522       include 'COMMON.LOCAL'
523       include 'COMMON.NAMES'
524       include 'COMMON.CHAIN'
525       include 'COMMON.FFIELD'
526       include 'COMMON.SBRIDGE'
527       include 'COMMON.HEADER'
528       include 'COMMON.CONTROL'
529       include 'COMMON.DBASE'
530       include 'COMMON.THREAD'
531       include 'COMMON.CONTACTS'
532       include 'COMMON.TORCNSTR'
533       include 'COMMON.TIME1'
534       include 'COMMON.BOUNDS'
535       include 'COMMON.MD'
536       include 'COMMON.SETUP'
537       character*4 sequence(maxres)
538       integer rescode
539       double precision x(maxvar)
540       character*256 pdbfile
541       character*320 weightcard
542       character*80 weightcard_t,ucase
543       dimension itype_pdb(maxres)
544       common /pizda/ itype_pdb
545       logical seq_comp,fail
546       double precision energia(0:n_ene)
547       integer ilen
548       external ilen
549 C
550 C Body
551 C
552 C Read weights of the subsequent energy terms.
553       call card_concat(weightcard)
554       call reada(weightcard,'WLONG',wlong,1.0D0)
555       call reada(weightcard,'WSC',wsc,wlong)
556       call reada(weightcard,'WSCP',wscp,wlong)
557       call reada(weightcard,'WELEC',welec,1.0D0)
558       call reada(weightcard,'WVDWPP',wvdwpp,welec)
559       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
560       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
561       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
562       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
563       call reada(weightcard,'WTURN3',wturn3,1.0D0)
564       call reada(weightcard,'WTURN4',wturn4,1.0D0)
565       call reada(weightcard,'WTURN6',wturn6,1.0D0)
566       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
567       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
568       call reada(weightcard,'WBOND',wbond,1.0D0)
569       call reada(weightcard,'WTOR',wtor,1.0D0)
570       call reada(weightcard,'WTORD',wtor_d,1.0D0)
571       call reada(weightcard,'WANG',wang,1.0D0)
572       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
573       call reada(weightcard,'SCAL14',scal14,0.4D0)
574       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
575       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
576       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
577       call reada(weightcard,'TEMP0',temp0,300.0d0)
578       if (index(weightcard,'SOFT').gt.0) ipot=6
579 C 12/1/95 Added weight for the multi-body term WCORR
580       call reada(weightcard,'WCORRH',wcorr,1.0D0)
581       if (wcorr4.gt.0.0d0) wcorr=wcorr4
582       weights(1)=wsc
583       weights(2)=wscp
584       weights(3)=welec
585       weights(4)=wcorr
586       weights(5)=wcorr5
587       weights(6)=wcorr6
588       weights(7)=wel_loc
589       weights(8)=wturn3
590       weights(9)=wturn4
591       weights(10)=wturn6
592       weights(11)=wang
593       weights(12)=wscloc
594       weights(13)=wtor
595       weights(14)=wtor_d
596       weights(15)=wstrain
597       weights(16)=wvdwpp
598       weights(17)=wbond
599       weights(18)=scal14
600       weights(21)=wsccor
601       if(me.eq.king.or..not.out1file)
602      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
603      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
604      &  wturn4,wturn6
605    10 format (/'Energy-term weights (unscaled):'//
606      & 'WSCC=   ',f10.6,' (SC-SC)'/
607      & 'WSCP=   ',f10.6,' (SC-p)'/
608      & 'WELEC=  ',f10.6,' (p-p electr)'/
609      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
610      & 'WBOND=  ',f10.6,' (stretching)'/
611      & 'WANG=   ',f10.6,' (bending)'/
612      & 'WSCLOC= ',f10.6,' (SC local)'/
613      & 'WTOR=   ',f10.6,' (torsional)'/
614      & 'WTORD=  ',f10.6,' (double torsional)'/
615      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
616      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
617      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
618      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
619      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
620      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
621      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
622      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
623      & 'WTURN6= ',f10.6,' (turns, 6th order)')
624       if(me.eq.king.or..not.out1file)then
625        if (wcorr4.gt.0.0d0) then
626         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
627      &   'between contact pairs of peptide groups'
628         write (iout,'(2(a,f5.3/))') 
629      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
630      &  'Range of quenching the correlation terms:',2*delt_corr 
631        else if (wcorr.gt.0.0d0) then
632         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
633      &   'between contact pairs of peptide groups'
634        endif
635        write (iout,'(a,f8.3)') 
636      &  'Scaling factor of 1,4 SC-p interactions:',scal14
637        write (iout,'(a,f8.3)') 
638      &  'General scaling factor of SC-p interactions:',scalscp
639       endif
640       r0_corr=cutoff_corr-delt_corr
641       do i=1,20
642         aad(i,1)=scalscp*aad(i,1)
643         aad(i,2)=scalscp*aad(i,2)
644         bad(i,1)=scalscp*bad(i,1)
645         bad(i,2)=scalscp*bad(i,2)
646       enddo
647       call rescale_weights(t_bath)
648       if(me.eq.king.or..not.out1file)
649      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
650      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
651      &  wturn4,wturn6
652    22 format (/'Energy-term weights (scaled):'//
653      & 'WSCC=   ',f10.6,' (SC-SC)'/
654      & 'WSCP=   ',f10.6,' (SC-p)'/
655      & 'WELEC=  ',f10.6,' (p-p electr)'/
656      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
657      & 'WBOND=  ',f10.6,' (stretching)'/
658      & 'WANG=   ',f10.6,' (bending)'/
659      & 'WSCLOC= ',f10.6,' (SC local)'/
660      & 'WTOR=   ',f10.6,' (torsional)'/
661      & 'WTORD=  ',f10.6,' (double torsional)'/
662      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
663      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
664      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
665      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
666      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
667      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
668      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
669      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
670      & 'WTURN6= ',f10.6,' (turns, 6th order)')
671       if(me.eq.king.or..not.out1file)
672      & write (iout,*) "Reference temperature for weights calculation:",
673      &  temp0
674       call reada(weightcard,"D0CM",d0cm,3.78d0)
675       call reada(weightcard,"AKCM",akcm,15.1d0)
676       call reada(weightcard,"AKTH",akth,11.0d0)
677       call reada(weightcard,"AKCT",akct,12.0d0)
678       call reada(weightcard,"V1SS",v1ss,-1.08d0)
679       call reada(weightcard,"V2SS",v2ss,7.61d0)
680       call reada(weightcard,"V3SS",v3ss,13.7d0)
681       call reada(weightcard,"EBR",ebr,-5.50D0)
682       if(me.eq.king.or..not.out1file) then
683        write (iout,*) "Parameters of the SS-bond potential:"
684        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
685      & " AKCT",akct
686        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
687        write (iout,*) "EBR",ebr
688 c       print *,'indpdb=',indpdb,' pdbref=',pdbref
689       endif
690       if (indpdb.gt.0 .or. pdbref) then
691         read(inp,'(a)') pdbfile
692         if(me.eq.king.or..not.out1file)
693      &   write (iout,'(2a)') 'PDB data will be read from file ',
694      &   pdbfile(:ilen(pdbfile))
695         open(ipdbin,file=pdbfile,status='old',err=33)
696         goto 34 
697   33    write (iout,'(a)') 'Error opening PDB file.'
698         stop
699   34    continue
700 c        print *,'Begin reading pdb data'
701         call readpdb
702 c        print *,'Finished reading pdb data'
703         if(me.eq.king.or..not.out1file)
704      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
705      &   ' nstart_sup=',nstart_sup
706         do i=1,nres
707           itype_pdb(i)=itype(i)
708         enddo
709         close (ipdbin)
710         nnt=nstart_sup
711         nct=nstart_sup+nsup-1
712         call contact(.false.,ncont_ref,icont_ref,co)
713
714         if (sideadd) then 
715          if(me.eq.king.or..not.out1file)
716      &    write(iout,*)'Adding sidechains'
717          maxsi=1000
718          do i=2,nres-1
719           iti=itype(i)
720           if (iti.ne.10 .and. itype(i).ne.21) then
721             nsi=0
722             fail=.true.
723             do while (fail.and.nsi.le.maxsi)
724               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
725               nsi=nsi+1
726             enddo
727             if(fail) write(iout,*)'Adding sidechain failed for res ',
728      &              i,' after ',nsi,' trials'
729           endif
730          enddo
731         endif  
732       endif
733       if (indpdb.eq.0) then
734 C Read sequence if not taken from the pdb file.
735         read (inp,*) nres
736 c        print *,'nres=',nres
737         if (iscode.gt.0) then
738           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
739         else
740           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
741         endif
742 C Convert sequence to numeric code
743         do i=1,nres
744           itype(i)=rescode(i,sequence(i),iscode)
745         enddo
746 C Assign initial virtual bond lengths
747         do i=2,nres
748           vbld(i)=vbl
749           vbld_inv(i)=vblinv
750         enddo
751         do i=2,nres-1
752           vbld(i+nres)=dsc(itype(i))
753           vbld_inv(i+nres)=dsc_inv(itype(i))
754 c          write (iout,*) "i",i," itype",itype(i),
755 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
756         enddo
757       endif 
758 c      print *,nres
759 c      print '(20i4)',(itype(i),i=1,nres)
760       do i=1,nres
761 #ifdef PROCOR
762         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
763 #else
764         if (itype(i).eq.21) then
765 #endif
766           itel(i)=0
767 #ifdef PROCOR
768         else if (itype(i+1).ne.20) then
769 #else
770         else if (itype(i).ne.20) then
771 #endif
772           itel(i)=1
773         else
774           itel(i)=2
775         endif  
776       enddo
777       if(me.eq.king.or..not.out1file)then
778        write (iout,*) "ITEL"
779        do i=1,nres-1
780          write (iout,*) i,itype(i),itel(i)
781        enddo
782        print *,'Call Read_Bridge.'
783       endif
784       call read_bridge
785 C 8/13/98 Set limits to generating the dihedral angles
786       do i=1,nres
787         phibound(1,i)=-pi
788         phibound(2,i)=pi
789       enddo
790       read (inp,*) ndih_constr
791       if (ndih_constr.gt.0) then
792         read (inp,*) ftors
793         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
794         if(me.eq.king.or..not.out1file)then
795          write (iout,*) 
796      &   'There are',ndih_constr,' constraints on phi angles.'
797          do i=1,ndih_constr
798           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
799          enddo
800         endif
801         do i=1,ndih_constr
802           phi0(i)=deg2rad*phi0(i)
803           drange(i)=deg2rad*drange(i)
804         enddo
805         if(me.eq.king.or..not.out1file)
806      &   write (iout,*) 'FTORS',ftors
807         do i=1,ndih_constr
808           ii = idih_constr(i)
809           phibound(1,ii) = phi0(i)-drange(i)
810           phibound(2,ii) = phi0(i)+drange(i)
811         enddo 
812       endif
813       nnt=1
814 #ifdef MPI
815       if (me.eq.king) then
816 #endif
817        write (iout,'(a)') 'Boundaries in phi angle sampling:'
818        do i=1,nres
819          write (iout,'(a3,i5,2f10.1)') 
820      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
821        enddo
822 #ifdef MP
823       endif
824 #endif
825       nct=nres
826 cd      print *,'NNT=',NNT,' NCT=',NCT
827       if (itype(1).eq.21) nnt=2
828       if (itype(nres).eq.21) nct=nct-1
829       if (pdbref) then
830         if(me.eq.king.or..not.out1file)
831      &   write (iout,'(a,i3)') 'nsup=',nsup
832         nstart_seq=nnt
833         if (nsup.le.(nct-nnt+1)) then
834           do i=0,nct-nnt+1-nsup
835             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
836               nstart_seq=nnt+i
837               goto 111
838             endif
839           enddo
840           write (iout,'(a)') 
841      &            'Error - sequences to be superposed do not match.'
842           stop
843         else
844           do i=0,nsup-(nct-nnt+1)
845             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
846      &      then
847               nstart_sup=nstart_sup+i
848               nsup=nct-nnt+1
849               goto 111
850             endif
851           enddo 
852           write (iout,'(a)') 
853      &            'Error - sequences to be superposed do not match.'
854         endif
855   111   continue
856         if (nsup.eq.0) nsup=nct-nnt
857         if (nstart_sup.eq.0) nstart_sup=nnt
858         if (nstart_seq.eq.0) nstart_seq=nnt
859         if(me.eq.king.or..not.out1file)  
860      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
861      &                 ' nstart_seq=',nstart_seq
862       endif
863 c--- Zscore rms -------
864       if (nz_start.eq.0) nz_start=nnt
865       if (nz_end.eq.0 .and. nsup.gt.0) then
866         nz_end=nnt+nsup-1
867       else if (nz_end.eq.0) then
868         nz_end=nct
869       endif
870       if(me.eq.king.or..not.out1file)then
871        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
872        write (iout,*) 'IZ_SC=',iz_sc
873       endif
874 c----------------------
875       call init_int_table
876       if (refstr) then
877         if (.not.pdbref) then
878           call read_angles(inp,*38)
879           goto 39
880    38     write (iout,'(a)') 'Error reading reference structure.'
881 #ifdef MPI
882           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
883           stop 'Error reading reference structure'
884 #endif
885    39     call chainbuild
886           call setup_var
887 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
888           nstart_sup=nnt
889           nstart_seq=nnt
890           nsup=nct-nnt+1
891           kkk=1
892           do i=1,2*nres
893             do j=1,3
894               cref(j,i,kkk)=c(j,i)
895             enddo
896           enddo
897           call contact(.true.,ncont_ref,icont_ref,co)
898         endif
899 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
900         call flush(iout)
901         if (constr_dist.gt.0) call read_dist_constr
902         write (iout,*) "After read_dist_constr nhpb",nhpb
903         call hpb_partition
904         if(me.eq.king.or..not.out1file)
905      &   write (iout,*) 'Contact order:',co
906         if (pdbref) then
907         if(me.eq.king.or..not.out1file)
908      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
909         do i=1,ncont_ref
910           do j=1,2
911             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
912           enddo
913           if(me.eq.king.or..not.out1file)
914      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
915      &     icont_ref(1,i),' ',
916      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
917         enddo
918         endif
919       endif
920       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
921      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
922      &    modecalc.ne.10) then
923 C If input structure hasn't been supplied from the PDB file read or generate
924 C initial geometry.
925         if (iranconf.eq.0 .and. .not. extconf) then
926           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
927      &     write (iout,'(a)') 'Initial geometry will be read in.'
928           if (read_cart) then
929             read(inp,'(8f10.5)',end=36,err=36)
930      &       ((c(l,k),l=1,3),k=1,nres),
931      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
932             write (iout,*) "Exit READ_CART"
933             write (iout,'(8f10.5)') 
934      &       ((c(l,k),l=1,3),k=1,nres),
935      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
936             call int_from_cart1(.true.)
937             write (iout,*) "Finish INT_TO_CART"
938             do i=1,nres-1
939               do j=1,3
940                 dc(j,i)=c(j,i+1)-c(j,i)
941                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
942               enddo
943             enddo
944             do i=nnt,nct
945               if (itype(i).ne.10 .and. itype(i).ne.21) then
946                 do j=1,3
947                   dc(j,i+nres)=c(j,i+nres)-c(j,i) 
948                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
949                 enddo
950               endif
951             enddo
952             return
953           else
954             call read_angles(inp,*36)
955           endif
956           goto 37
957    36     write (iout,'(a)') 'Error reading angle file.'
958 #ifdef MPI
959           call mpi_finalize( MPI_COMM_WORLD,IERR )
960 #endif
961           stop 'Error reading angle file.'
962    37     continue 
963         else if (extconf) then
964          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
965      &    write (iout,'(a)') 'Extended chain initial geometry.'
966          do i=3,nres
967           theta(i)=90d0*deg2rad
968          enddo
969          do i=4,nres
970           phi(i)=180d0*deg2rad
971          enddo
972          do i=2,nres-1
973           alph(i)=110d0*deg2rad
974          enddo
975          do i=2,nres-1
976           omeg(i)=-120d0*deg2rad
977          enddo
978         else
979           if(me.eq.king.or..not.out1file)
980      &     write (iout,'(a)') 'Random-generated initial geometry.'
981
982
983 #ifdef MPI
984           if (me.eq.king  .or. fg_rank.eq.0 .and. (
985      &           modecalc.eq.12 .or. modecalc.eq.14) ) then  
986 #endif
987             do itrial=1,100
988               itmp=1
989               call gen_rand_conf(itmp,*30)
990               goto 40
991    30         write (iout,*) 'Failed to generate random conformation',
992      &          ', itrial=',itrial
993               write (*,*) 'Processor:',me,
994      &          ' Failed to generate random conformation',
995      &          ' itrial=',itrial
996               call intout
997
998 #ifdef AIX
999               call flush_(iout)
1000 #else
1001               call flush(iout)
1002 #endif
1003             enddo
1004             write (iout,'(a,i3,a)') 'Processor:',me,
1005      &        ' error in generating random conformation.'
1006             write (*,'(a,i3,a)') 'Processor:',me,
1007      &        ' error in generating random conformation.'
1008             call flush(iout)
1009 #ifdef MPI
1010             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)            
1011    40       continue
1012           endif
1013 #else
1014           do itrial=1,100
1015             itmp=1
1016             call gen_rand_conf(itmp,*30)
1017             goto 40
1018    30       write (iout,*) 'Failed to generate random conformation',
1019      &        ', itrial=',itrial
1020             write (*,*) 'Failed to generate random conformation',
1021      &        ', itrial=',itrial
1022           enddo
1023           write (iout,'(a,i3,a)') 'Processor:',me,
1024      &      ' error in generating random conformation.'
1025           write (*,'(a,i3,a)') 'Processor:',me,
1026      &      ' error in generating random conformation.'
1027           stop
1028    40     continue
1029 #endif
1030         endif
1031       elseif (modecalc.eq.4) then
1032         read (inp,'(a)') intinname
1033         open (intin,file=intinname,status='old',err=333)
1034         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1035      &  write (iout,'(a)') 'intinname',intinname
1036         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1037         goto 334
1038   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1039 #ifdef MPI 
1040         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1041 #endif   
1042         stop 'Error opening angle file.' 
1043   334   continue
1044
1045       endif 
1046 C Generate distance constraints, if the PDB structure is to be regularized. 
1047       if (nthread.gt.0) then
1048         call read_threadbase
1049       endif
1050       call setup_var
1051       if (me.eq.king .or. .not. out1file)
1052      & call intout
1053       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1054         write (iout,'(/a,i3,a)') 
1055      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1056         write (iout,'(20i4)') (iss(i),i=1,ns)
1057         write (iout,'(/a/)') 'Pre-formed links are:' 
1058         do i=1,nss
1059           i1=ihpb(i)-nres
1060           i2=jhpb(i)-nres
1061           it1=itype(i1)
1062           it2=itype(i2)
1063           if (me.eq.king.or..not.out1file)
1064      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1065      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1066      &    ebr,forcon(i)
1067         enddo
1068         write (iout,'(a)')
1069       endif
1070       if (i2ndstr.gt.0) call secstrp2dihc
1071 c      call geom_to_var(nvar,x)
1072 c      call etotal(energia(0))
1073 c      call enerprint(energia(0))
1074 c      call briefout(0,etot)
1075 c      stop
1076 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1077 cd    write (iout,'(a)') 'Variable list:'
1078 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1079 #ifdef MPI
1080       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1081      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
1082      &  'Processor',myrank,': end reading molecular data.'
1083 #endif
1084       return
1085       end
1086 c--------------------------------------------------------------------------
1087       logical function seq_comp(itypea,itypeb,length)
1088       implicit none
1089       integer length,itypea(length),itypeb(length)
1090       integer i
1091       do i=1,length
1092         if (itypea(i).ne.itypeb(i)) then
1093           seq_comp=.false.
1094           return
1095         endif
1096       enddo
1097       seq_comp=.true.
1098       return
1099       end
1100 c-----------------------------------------------------------------------------
1101       subroutine read_bridge
1102 C Read information about disulfide bridges.
1103       implicit real*8 (a-h,o-z)
1104       include 'DIMENSIONS'
1105 #ifdef MPI
1106       include 'mpif.h'
1107 #endif
1108       include 'COMMON.IOUNITS'
1109       include 'COMMON.GEO'
1110       include 'COMMON.VAR'
1111       include 'COMMON.INTERACT'
1112       include 'COMMON.LOCAL'
1113       include 'COMMON.NAMES'
1114       include 'COMMON.CHAIN'
1115       include 'COMMON.FFIELD'
1116       include 'COMMON.SBRIDGE'
1117       include 'COMMON.HEADER'
1118       include 'COMMON.CONTROL'
1119       include 'COMMON.DBASE'
1120       include 'COMMON.THREAD'
1121       include 'COMMON.TIME1'
1122       include 'COMMON.SETUP'
1123 C Read bridging residues.
1124       read (inp,*) ns,(iss(i),i=1,ns)
1125 c      print *,'ns=',ns
1126       if(me.eq.king.or..not.out1file)
1127      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1128 C Check whether the specified bridging residues are cystines.
1129       do i=1,ns
1130         if (itype(iss(i)).ne.1) then
1131           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') 
1132      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1133      &   ' can form a disulfide bridge?!!!'
1134           write (*,'(2a,i3,a)') 
1135      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1136      &   ' can form a disulfide bridge?!!!'
1137 #ifdef MPI
1138          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1139          stop
1140 #endif
1141         endif
1142       enddo
1143 C Read preformed bridges.
1144       if (ns.gt.0) then
1145       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1146       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1147       if (nss.gt.0) then
1148         nhpb=nss
1149 C Check if the residues involved in bridges are in the specified list of
1150 C bridging residues.
1151         do i=1,nss
1152           do j=1,i-1
1153             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1154      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1155               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1156      &      ' contains residues present in other pairs.'
1157               write (*,'(a,i3,a)') 'Disulfide pair',i,
1158      &      ' contains residues present in other pairs.'
1159 #ifdef MPI
1160               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1161               stop 
1162 #endif
1163             endif
1164           enddo
1165           do j=1,ns
1166             if (ihpb(i).eq.iss(j)) goto 10
1167           enddo
1168           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1169    10     continue
1170           do j=1,ns
1171             if (jhpb(i).eq.iss(j)) goto 20
1172           enddo
1173           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1174    20     continue
1175           dhpb(i)=dbr
1176           forcon(i)=fbr
1177         enddo
1178         do i=1,nss
1179           ihpb(i)=ihpb(i)+nres
1180           jhpb(i)=jhpb(i)+nres
1181         enddo
1182       endif
1183       endif
1184       return
1185       end
1186 c----------------------------------------------------------------------------
1187       subroutine read_x(kanal,*)
1188       implicit real*8 (a-h,o-z)
1189       include 'DIMENSIONS'
1190       include 'COMMON.GEO'
1191       include 'COMMON.VAR'
1192       include 'COMMON.CHAIN'
1193       include 'COMMON.IOUNITS'
1194       include 'COMMON.CONTROL'
1195       include 'COMMON.LOCAL'
1196       include 'COMMON.INTERACT'
1197 c Read coordinates from input
1198 c
1199       read(kanal,'(8f10.5)',end=10,err=10)
1200      &  ((c(l,k),l=1,3),k=1,nres),
1201      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1202       do j=1,3
1203         c(j,nres+1)=c(j,1)
1204         c(j,2*nres)=c(j,nres)
1205       enddo
1206       call int_from_cart1(.false.)
1207       do i=1,nres-1
1208         do j=1,3
1209           dc(j,i)=c(j,i+1)-c(j,i)
1210           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1211         enddo
1212       enddo
1213       do i=nnt,nct
1214         if (itype(i).ne.10 .and. itype(i).ne.21) then
1215           do j=1,3
1216             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1217             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1218           enddo
1219         endif
1220       enddo
1221
1222       return
1223    10 return1
1224       end
1225 c----------------------------------------------------------------------------
1226       subroutine read_threadbase
1227       implicit real*8 (a-h,o-z)
1228       include 'DIMENSIONS'
1229       include 'COMMON.IOUNITS'
1230       include 'COMMON.GEO'
1231       include 'COMMON.VAR'
1232       include 'COMMON.INTERACT'
1233       include 'COMMON.LOCAL'
1234       include 'COMMON.NAMES'
1235       include 'COMMON.CHAIN'
1236       include 'COMMON.FFIELD'
1237       include 'COMMON.SBRIDGE'
1238       include 'COMMON.HEADER'
1239       include 'COMMON.CONTROL'
1240       include 'COMMON.DBASE'
1241       include 'COMMON.THREAD'
1242       include 'COMMON.TIME1'
1243 C Read pattern database for threading.
1244       read (icbase,*) nseq
1245       do i=1,nseq
1246         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1247      &   nres_base(2,i),nres_base(3,i)
1248         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1249      &   nres_base(1,i))
1250 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1251 c    &   nres_base(2,i),nres_base(3,i)
1252 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1253 c    &   nres_base(1,i))
1254       enddo
1255       close (icbase)
1256       if (weidis.eq.0.0D0) weidis=0.1D0
1257       do i=nnt,nct
1258         do j=i+2,nct
1259           nhpb=nhpb+1
1260           ihpb(nhpb)=i
1261           jhpb(nhpb)=j
1262           forcon(nhpb)=weidis
1263         enddo
1264       enddo 
1265       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1266       write (iout,'(a,i5)') 'nexcl: ',nexcl
1267       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1268       return
1269       end
1270 c------------------------------------------------------------------------------
1271       subroutine setup_var
1272       implicit real*8 (a-h,o-z)
1273       include 'DIMENSIONS'
1274       include 'COMMON.IOUNITS'
1275       include 'COMMON.GEO'
1276       include 'COMMON.VAR'
1277       include 'COMMON.INTERACT'
1278       include 'COMMON.LOCAL'
1279       include 'COMMON.NAMES'
1280       include 'COMMON.CHAIN'
1281       include 'COMMON.FFIELD'
1282       include 'COMMON.SBRIDGE'
1283       include 'COMMON.HEADER'
1284       include 'COMMON.CONTROL'
1285       include 'COMMON.DBASE'
1286       include 'COMMON.THREAD'
1287       include 'COMMON.TIME1'
1288 C Set up variable list.
1289       ntheta=nres-2
1290       nphi=nres-3
1291       nvar=ntheta+nphi
1292       nside=0
1293       do i=2,nres-1
1294         if (itype(i).ne.10 .and. itype(i).ne.21) then
1295           nside=nside+1
1296           ialph(i,1)=nvar+nside
1297           ialph(nside,2)=i
1298         endif
1299       enddo
1300       if (indphi.gt.0) then
1301         nvar=nphi
1302       else if (indback.gt.0) then
1303         nvar=nphi+ntheta
1304       else
1305         nvar=nvar+2*nside
1306       endif
1307 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1308       return
1309       end
1310 c----------------------------------------------------------------------------
1311       subroutine gen_dist_constr
1312 C Generate CA distance constraints.
1313       implicit real*8 (a-h,o-z)
1314       include 'DIMENSIONS'
1315       include 'COMMON.IOUNITS'
1316       include 'COMMON.GEO'
1317       include 'COMMON.VAR'
1318       include 'COMMON.INTERACT'
1319       include 'COMMON.LOCAL'
1320       include 'COMMON.NAMES'
1321       include 'COMMON.CHAIN'
1322       include 'COMMON.FFIELD'
1323       include 'COMMON.SBRIDGE'
1324       include 'COMMON.HEADER'
1325       include 'COMMON.CONTROL'
1326       include 'COMMON.DBASE'
1327       include 'COMMON.THREAD'
1328       include 'COMMON.TIME1'
1329       dimension itype_pdb(maxres)
1330       common /pizda/ itype_pdb
1331       character*2 iden
1332 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1333 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1334 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1335 cd     & ' nsup',nsup
1336       do i=nstart_sup,nstart_sup+nsup-1
1337 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1338 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1339         do j=i+2,nstart_sup+nsup-1
1340           nhpb=nhpb+1
1341           ihpb(nhpb)=i+nstart_seq-nstart_sup
1342           jhpb(nhpb)=j+nstart_seq-nstart_sup
1343           forcon(nhpb)=weidis
1344           dhpb(nhpb)=dist(i,j)
1345         enddo
1346       enddo 
1347 cd      write (iout,'(a)') 'Distance constraints:' 
1348 cd      do i=nss+1,nhpb
1349 cd        ii=ihpb(i)
1350 cd        jj=jhpb(i)
1351 cd        iden='CA'
1352 cd        if (ii.gt.nres) then
1353 cd          iden='SC'
1354 cd          ii=ii-nres
1355 cd          jj=jj-nres
1356 cd        endif
1357 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') 
1358 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1359 cd     &  dhpb(i),forcon(i)
1360 cd      enddo
1361       return
1362       end
1363 c----------------------------------------------------------------------------
1364       subroutine map_read
1365       implicit real*8 (a-h,o-z)
1366       include 'DIMENSIONS'
1367       include 'COMMON.MAP'
1368       include 'COMMON.IOUNITS'
1369       character*3 angid(4) /'THE','PHI','ALP','OME'/
1370       character*80 mapcard,ucase
1371       do imap=1,nmap
1372         read (inp,'(a)') mapcard
1373         mapcard=ucase(mapcard)
1374         if (index(mapcard,'PHI').gt.0) then
1375           kang(imap)=1
1376         else if (index(mapcard,'THE').gt.0) then
1377           kang(imap)=2
1378         else if (index(mapcard,'ALP').gt.0) then
1379           kang(imap)=3
1380         else if (index(mapcard,'OME').gt.0) then
1381           kang(imap)=4
1382         else
1383           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1384           stop 'Error - illegal variable spec in MAP card.'
1385         endif
1386         call readi (mapcard,'RES1',res1(imap),0)
1387         call readi (mapcard,'RES2',res2(imap),0)
1388         if (res1(imap).eq.0) then
1389           res1(imap)=res2(imap)
1390         else if (res2(imap).eq.0) then
1391           res2(imap)=res1(imap)
1392         endif
1393         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1394           write (iout,'(a)') 
1395      &    'Error - illegal definition of variable group in MAP.'
1396           stop 'Error - illegal definition of variable group in MAP.'
1397         endif
1398         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1399         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1400         call readi(mapcard,'NSTEP',nstep(imap),0)
1401         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1402           write (iout,'(a)') 
1403      &     'Illegal boundary and/or step size specification in MAP.'
1404           stop 'Illegal boundary and/or step size specification in MAP.'
1405         endif
1406       enddo ! imap
1407       return
1408       end 
1409 c----------------------------------------------------------------------------
1410       subroutine csaread
1411       implicit real*8 (a-h,o-z)
1412       include 'DIMENSIONS'
1413       include 'COMMON.IOUNITS'
1414       include 'COMMON.GEO'
1415       include 'COMMON.CSA'
1416       include 'COMMON.BANK'
1417       include 'COMMON.CONTROL'
1418       character*80 ucase
1419       character*620 mcmcard
1420       call card_concat(mcmcard)
1421
1422       call readi(mcmcard,'NCONF',nconf,50)
1423       call readi(mcmcard,'NADD',nadd,0)
1424       call readi(mcmcard,'JSTART',jstart,1)
1425       call readi(mcmcard,'JEND',jend,1)
1426       call readi(mcmcard,'NSTMAX',nstmax,500000)
1427       call readi(mcmcard,'N0',n0,1)
1428       call readi(mcmcard,'N1',n1,6)
1429       call readi(mcmcard,'N2',n2,4)
1430       call readi(mcmcard,'N3',n3,0)
1431       call readi(mcmcard,'N4',n4,0)
1432       call readi(mcmcard,'N5',n5,0)
1433       call readi(mcmcard,'N6',n6,10)
1434       call readi(mcmcard,'N7',n7,0)
1435       call readi(mcmcard,'N8',n8,0)
1436       call readi(mcmcard,'N9',n9,0)
1437       call readi(mcmcard,'N14',n14,0)
1438       call readi(mcmcard,'N15',n15,0)
1439       call readi(mcmcard,'N16',n16,0)
1440       call readi(mcmcard,'N17',n17,0)
1441       call readi(mcmcard,'N18',n18,0)
1442
1443       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1444
1445       call readi(mcmcard,'NDIFF',ndiff,2)
1446       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1447       call readi(mcmcard,'IS1',is1,1)
1448       call readi(mcmcard,'IS2',is2,8)
1449       call readi(mcmcard,'NRAN0',nran0,4)
1450       call readi(mcmcard,'NRAN1',nran1,2)
1451       call readi(mcmcard,'IRR',irr,1)
1452       call readi(mcmcard,'NSEED',nseed,20)
1453       call readi(mcmcard,'NTOTAL',ntotal,10000)
1454       call reada(mcmcard,'CUT1',cut1,2.0d0)
1455       call reada(mcmcard,'CUT2',cut2,5.0d0)
1456       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1457       call readi(mcmcard,'ICMAX',icmax,3)
1458       call readi(mcmcard,'IRESTART',irestart,0)
1459 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1460       ntbankm=0
1461 c!bankt
1462       call reada(mcmcard,'DELE',dele,20.0d0)
1463       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1464       call readi(mcmcard,'IREF',iref,0)
1465       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1466       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1467       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1468       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1469       write (iout,*) "NCONF_IN",nconf_in
1470       return
1471       end
1472 c----------------------------------------------------------------------------
1473 cfmc      subroutine mcmfread
1474 cfmc      implicit real*8 (a-h,o-z)
1475 cfmc      include 'DIMENSIONS'
1476 cfmc      include 'COMMON.MCMF'
1477 cfmc      include 'COMMON.IOUNITS'
1478 cfmc      include 'COMMON.GEO'
1479 cfmc      character*80 ucase
1480 cfmc      character*620 mcmcard
1481 cfmc      call card_concat(mcmcard)
1482 cfmc
1483 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1484 cfmc      write(iout,*)'MAXRANT=',maxrant
1485 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1486 cfmc      write(iout,*)'MAXFAM=',maxfam
1487 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1488 cfmc      write(iout,*)'NNET1=',nnet1
1489 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1490 cfmc      write(iout,*)'NNET2=',nnet2
1491 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1492 cfmc      write(iout,*)'NNET3=',nnet3
1493 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1494 cfmc      write(iout,*)'ILASTT=',ilastt
1495 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1496 cfmc      write(iout,*)'MAXSTR=',maxstr
1497 cfmc      maxstr_f=maxstr/maxfam
1498 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1499 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1500 cfmc      write(iout,*)'NMCMF=',nmcmf
1501 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1502 cfmc      write(iout,*)'IFOCUS=',ifocus
1503 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1504 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1505 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1506 cfmc      write(iout,*)'INTPRT=',intprt
1507 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1508 cfmc      write(iout,*)'IPRT=',iprt
1509 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1510 cfmc      write(iout,*)'IMAXTR=',imaxtr
1511 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1512 cfmc      write(iout,*)'MAXEVEN=',maxeven
1513 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1514 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1515 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1516 cfmc      write(iout,*)'INIMIN=',inimin
1517 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1518 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1519 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1520 cfmc      write(iout,*)'NTHREAD=',nthread
1521 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1522 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1523 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1524 cfmc      write(iout,*)'MAXPERT=',maxpert
1525 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1526 cfmc      write(iout,*)'IRMSD=',irmsd
1527 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1528 cfmc      write(iout,*)'DENEMIN=',denemin
1529 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1530 cfmc      write(iout,*)'RCUT1S=',rcut1s
1531 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1532 cfmc      write(iout,*)'RCUT1E=',rcut1e
1533 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1534 cfmc      write(iout,*)'RCUT2S=',rcut2s
1535 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1536 cfmc      write(iout,*)'RCUT2E=',rcut2e
1537 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1538 cfmc      write(iout,*)'DPERT1=',d_pert1
1539 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1540 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1541 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1542 cfmc      write(iout,*)'DPERT2=',d_pert2
1543 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1544 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1545 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1546 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1547 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1548 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1549 cfmc      d_pert1=deg2rad*d_pert1
1550 cfmc      d_pert1a=deg2rad*d_pert1a
1551 cfmc      d_pert2=deg2rad*d_pert2
1552 cfmc      d_pert2a=deg2rad*d_pert2a
1553 cfmc      d_pert2b=deg2rad*d_pert2b
1554 cfmc      d_pert2c=deg2rad*d_pert2c
1555 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1556 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1557 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1558 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1559 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1560 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1561 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1562 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1563 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1564 cfmc      write(iout,*)'RCUTINI=',rcutini
1565 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1566 cfmc      write(iout,*)'GRAT=',grat
1567 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1568 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1569 cfmc
1570 cfmc      return
1571 cfmc      end 
1572 c----------------------------------------------------------------------------
1573       subroutine mcmread
1574       implicit real*8 (a-h,o-z)
1575       include 'DIMENSIONS'
1576       include 'COMMON.MCM'
1577       include 'COMMON.MCE'
1578       include 'COMMON.IOUNITS'
1579       character*80 ucase
1580       character*320 mcmcard
1581       call card_concat(mcmcard)
1582       call readi(mcmcard,'MAXACC',maxacc,100)
1583       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1584       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1585       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1586       call readi(mcmcard,'MAXREPM',maxrepm,200)
1587       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1588       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1589       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1590       call reada(mcmcard,'E_UP',e_up,5.0D0)
1591       call reada(mcmcard,'DELTE',delte,0.1D0)
1592       call readi(mcmcard,'NSWEEP',nsweep,5)
1593       call readi(mcmcard,'NSTEPH',nsteph,0)
1594       call readi(mcmcard,'NSTEPC',nstepc,0)
1595       call reada(mcmcard,'TMIN',tmin,298.0D0)
1596       call reada(mcmcard,'TMAX',tmax,298.0D0)
1597       call readi(mcmcard,'NWINDOW',nwindow,0)
1598       call readi(mcmcard,'PRINT_MC',print_mc,0)
1599       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1600       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1601       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1602       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1603       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1604       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1605       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1606       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1607       if (nwindow.gt.0) then
1608         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1609         do i=1,nwindow
1610           winlen(i)=winend(i)-winstart(i)+1
1611         enddo
1612       endif
1613       if (tmax.lt.tmin) tmax=tmin
1614       if (tmax.eq.tmin) then
1615         nstepc=0
1616         nsteph=0
1617       endif
1618       if (nstepc.gt.0 .and. nsteph.gt.0) then
1619         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0)) 
1620         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0)) 
1621       endif
1622 C Probabilities of different move types
1623       sumpro_type(0)=0.0D0
1624       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1625       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1626       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1627       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1628       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1629       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1630       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1631       do i=1,MaxMoveType
1632         print *,'i',i,' sumprotype',sumpro_type(i)
1633         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1634         print *,'i',i,' sumprotype',sumpro_type(i)
1635       enddo
1636       return
1637       end 
1638 c----------------------------------------------------------------------------
1639       subroutine read_minim
1640       implicit real*8 (a-h,o-z)
1641       include 'DIMENSIONS'
1642       include 'COMMON.MINIM'
1643       include 'COMMON.IOUNITS'
1644       character*80 ucase
1645       character*320 minimcard
1646       call card_concat(minimcard)
1647       call readi(minimcard,'MAXMIN',maxmin,2000)
1648       call readi(minimcard,'MAXFUN',maxfun,5000)
1649       call readi(minimcard,'MINMIN',minmin,maxmin)
1650       call readi(minimcard,'MINFUN',minfun,maxmin)
1651       call reada(minimcard,'TOLF',tolf,1.0D-2)
1652       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1653       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1654       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1655       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1656       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1657      &         'Options in energy minimization:'
1658       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1659      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1660      & 'MinMin:',MinMin,' MinFun:',MinFun,
1661      & ' TolF:',TolF,' RTolF:',RTolF
1662       return
1663       end
1664 c----------------------------------------------------------------------------
1665       subroutine read_angles(kanal,*)
1666       implicit real*8 (a-h,o-z)
1667       include 'DIMENSIONS'
1668       include 'COMMON.GEO'
1669       include 'COMMON.VAR'
1670       include 'COMMON.CHAIN'
1671       include 'COMMON.IOUNITS'
1672       include 'COMMON.CONTROL'
1673 c Read angles from input 
1674 c
1675        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1676        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1677        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1678        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1679
1680        do i=1,nres
1681 c 9/7/01 avoid 180 deg valence angle
1682         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1683 c
1684         theta(i)=deg2rad*theta(i)
1685         phi(i)=deg2rad*phi(i)
1686         alph(i)=deg2rad*alph(i)
1687         omeg(i)=deg2rad*omeg(i)
1688        enddo
1689       return
1690    10 return1
1691       end
1692 c----------------------------------------------------------------------------
1693       subroutine reada(rekord,lancuch,wartosc,default)
1694       implicit none
1695       character*(*) rekord,lancuch
1696       double precision wartosc,default
1697       integer ilen,iread
1698       external ilen
1699       iread=index(rekord,lancuch)
1700       if (iread.eq.0) then
1701         wartosc=default 
1702         return
1703       endif   
1704       iread=iread+ilen(lancuch)+1
1705       read (rekord(iread:),*,err=10,end=10) wartosc
1706       return
1707   10  wartosc=default
1708       return
1709       end
1710 c----------------------------------------------------------------------------
1711       subroutine readi(rekord,lancuch,wartosc,default)
1712       implicit none
1713       character*(*) rekord,lancuch
1714       integer wartosc,default
1715       integer ilen,iread
1716       external ilen
1717       iread=index(rekord,lancuch)
1718       if (iread.eq.0) then
1719         wartosc=default 
1720         return
1721       endif   
1722       iread=iread+ilen(lancuch)+1
1723       read (rekord(iread:),*,err=10,end=10) wartosc
1724       return
1725   10  wartosc=default
1726       return
1727       end
1728 c----------------------------------------------------------------------------
1729       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1730       implicit none
1731       integer dim,i
1732       integer tablica(dim),default
1733       character*(*) rekord,lancuch
1734       character*80 aux
1735       integer ilen,iread
1736       external ilen
1737       do i=1,dim
1738         tablica(i)=default
1739       enddo
1740       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1741       if (iread.eq.0) return
1742       iread=iread+ilen(lancuch)+1
1743       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1744    10 return
1745       end
1746 c----------------------------------------------------------------------------
1747       subroutine multreada(rekord,lancuch,tablica,dim,default)
1748       implicit none
1749       integer dim,i
1750       double precision tablica(dim),default
1751       character*(*) rekord,lancuch
1752       character*80 aux
1753       integer ilen,iread
1754       external ilen
1755       do i=1,dim
1756         tablica(i)=default
1757       enddo
1758       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1759       if (iread.eq.0) return
1760       iread=iread+ilen(lancuch)+1
1761       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1762    10 return
1763       end
1764 c----------------------------------------------------------------------------
1765       subroutine openunits
1766       implicit real*8 (a-h,o-z)
1767       include 'DIMENSIONS'    
1768 #ifdef MPI
1769       include 'mpif.h'
1770       character*16 form,nodename
1771       integer nodelen
1772 #endif
1773       include 'COMMON.SETUP'
1774       include 'COMMON.IOUNITS'
1775       include 'COMMON.MD'
1776       include 'COMMON.CONTROL'
1777       integer lenpre,lenpot,ilen,lentmp
1778       external ilen
1779       character*3 out1file_text,ucase
1780       character*3 ll
1781       external ucase
1782 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1783       call getenv_loc("PREFIX",prefix)
1784       pref_orig = prefix
1785       call getenv_loc("POT",pot)
1786       call getenv_loc("DIRTMP",tmpdir)
1787       call getenv_loc("CURDIR",curdir)
1788       call getenv_loc("OUT1FILE",out1file_text)
1789 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1790       out1file_text=ucase(out1file_text)
1791       if (out1file_text(1:1).eq."Y") then
1792         out1file=.true.
1793       else 
1794         out1file=fg_rank.gt.0
1795       endif
1796       lenpre=ilen(prefix)
1797       lenpot=ilen(pot)
1798       lentmp=ilen(tmpdir)
1799       if (lentmp.gt.0) then
1800           write (*,'(80(1h!))')
1801           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1802           write (*,'(80(1h!))')
1803           write (*,*)"All output files will be on node /tmp directory." 
1804 #ifdef MPI
1805         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1806         if (me.eq.king) then
1807           write (*,*) "The master node is ",nodename
1808         else if (fg_rank.eq.0) then
1809           write (*,*) "I am the CG slave node ",nodename
1810         else 
1811           write (*,*) "I am the FG slave node ",nodename
1812         endif
1813 #endif
1814         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1815         lenpre = lentmp+lenpre+1
1816       endif
1817       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1818 C Get the names and open the input files
1819 #if defined(WINIFL) || defined(WINPGI)
1820       open(1,file=pref_orig(:ilen(pref_orig))//
1821      &  '.inp',status='old',readonly,shared)
1822        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1823 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1824 C Get parameter filenames and open the parameter files.
1825       call getenv_loc('BONDPAR',bondname)
1826       open (ibond,file=bondname,status='old',readonly,shared)
1827       call getenv_loc('THETPAR',thetname)
1828       open (ithep,file=thetname,status='old',readonly,shared)
1829       call getenv_loc('ROTPAR',rotname)
1830       open (irotam,file=rotname,status='old',readonly,shared)
1831       call getenv_loc('TORPAR',torname)
1832       open (itorp,file=torname,status='old',readonly,shared)
1833       call getenv_loc('TORDPAR',tordname)
1834       open (itordp,file=tordname,status='old',readonly,shared)
1835       call getenv_loc('FOURIER',fouriername)
1836       open (ifourier,file=fouriername,status='old',readonly,shared)
1837       call getenv_loc('ELEPAR',elename)
1838       open (ielep,file=elename,status='old',readonly,shared)
1839       call getenv_loc('SIDEPAR',sidename)
1840       open (isidep,file=sidename,status='old',readonly,shared)
1841 #elif (defined CRAY) || (defined AIX)
1842       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1843      &  action='read')
1844 c      print *,"Processor",myrank," opened file 1" 
1845       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1846 c      print *,"Processor",myrank," opened file 9" 
1847 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1848 C Get parameter filenames and open the parameter files.
1849       call getenv_loc('BONDPAR',bondname)
1850       open (ibond,file=bondname,status='old',action='read')
1851 c      print *,"Processor",myrank," opened file IBOND" 
1852       call getenv_loc('THETPAR',thetname)
1853       open (ithep,file=thetname,status='old',action='read')
1854 c      print *,"Processor",myrank," opened file ITHEP" 
1855       call getenv_loc('ROTPAR',rotname)
1856       open (irotam,file=rotname,status='old',action='read')
1857 c      print *,"Processor",myrank," opened file IROTAM" 
1858       call getenv_loc('TORPAR',torname)
1859       open (itorp,file=torname,status='old',action='read')
1860 c      print *,"Processor",myrank," opened file ITORP" 
1861       call getenv_loc('TORDPAR',tordname)
1862       open (itordp,file=tordname,status='old',action='read')
1863 c      print *,"Processor",myrank," opened file ITORDP" 
1864       call getenv_loc('SCCORPAR',sccorname)
1865       open (isccor,file=sccorname,status='old',action='read')
1866 c      print *,"Processor",myrank," opened file ISCCOR" 
1867       call getenv_loc('FOURIER',fouriername)
1868       open (ifourier,file=fouriername,status='old',action='read')
1869 c      print *,"Processor",myrank," opened file IFOURIER" 
1870       call getenv_loc('ELEPAR',elename)
1871       open (ielep,file=elename,status='old',action='read')
1872 c      print *,"Processor",myrank," opened file IELEP" 
1873       call getenv_loc('SIDEPAR',sidename)
1874       open (isidep,file=sidename,status='old',action='read')
1875 c      print *,"Processor",myrank," opened file ISIDEP" 
1876 c      print *,"Processor",myrank," opened parameter files" 
1877 #elif (defined G77)
1878       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1879       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1880 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1881 C Get parameter filenames and open the parameter files.
1882       call getenv_loc('BONDPAR',bondname)
1883       open (ibond,file=bondname,status='old')
1884       call getenv_loc('THETPAR',thetname)
1885       open (ithep,file=thetname,status='old')
1886       call getenv_loc('ROTPAR',rotname)
1887       open (irotam,file=rotname,status='old')
1888       call getenv_loc('TORPAR',torname)
1889       open (itorp,file=torname,status='old')
1890       call getenv_loc('TORDPAR',tordname)
1891       open (itordp,file=tordname,status='old')
1892       call getenv_loc('SCCORPAR',sccorname)
1893       open (isccor,file=sccorname,status='old')
1894       call getenv_loc('FOURIER',fouriername)
1895       open (ifourier,file=fouriername,status='old')
1896       call getenv_loc('ELEPAR',elename)
1897       open (ielep,file=elename,status='old')
1898       call getenv_loc('SIDEPAR',sidename)
1899       open (isidep,file=sidename,status='old')
1900 #else
1901       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1902      &  readonly)
1903        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1904 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1905 C Get parameter filenames and open the parameter files.
1906       call getenv_loc('BONDPAR',bondname)
1907       open (ibond,file=bondname,status='old',readonly)
1908       call getenv_loc('THETPAR',thetname)
1909       open (ithep,file=thetname,status='old',readonly)
1910       call getenv_loc('ROTPAR',rotname)
1911       open (irotam,file=rotname,status='old',readonly)
1912       call getenv_loc('TORPAR',torname)
1913       open (itorp,file=torname,status='old',readonly)
1914       call getenv_loc('TORDPAR',tordname)
1915       open (itordp,file=tordname,status='old',readonly)
1916       call getenv_loc('SCCORPAR',sccorname)
1917       open (isccor,file=sccorname,status='old',readonly)
1918       call getenv_loc('FOURIER',fouriername)
1919       open (ifourier,file=fouriername,status='old',readonly)
1920       call getenv_loc('ELEPAR',elename)
1921       open (ielep,file=elename,status='old',readonly)
1922       call getenv_loc('SIDEPAR',sidename)
1923       open (isidep,file=sidename,status='old',readonly)
1924 #endif
1925 #ifndef OLDSCP
1926 C
1927 C 8/9/01 In the newest version SCp interaction constants are read from a file
1928 C Use -DOLDSCP to use hard-coded constants instead.
1929 C
1930       call getenv_loc('SCPPAR',scpname)
1931 #if defined(WINIFL) || defined(WINPGI)
1932       open (iscpp,file=scpname,status='old',readonly,shared)
1933 #elif (defined CRAY)  || (defined AIX)
1934       open (iscpp,file=scpname,status='old',action='read')
1935 #elif (defined G77)
1936       open (iscpp,file=scpname,status='old')
1937 #else
1938       open (iscpp,file=scpname,status='old',readonly)
1939 #endif
1940 #endif
1941       call getenv_loc('PATTERN',patname)
1942 #if defined(WINIFL) || defined(WINPGI)
1943       open (icbase,file=patname,status='old',readonly,shared)
1944 #elif (defined CRAY)  || (defined AIX)
1945       open (icbase,file=patname,status='old',action='read')
1946 #elif (defined G77)
1947       open (icbase,file=patname,status='old')
1948 #else
1949       open (icbase,file=patname,status='old',readonly)
1950 #endif
1951 #ifdef MPI
1952 C Open output file only for CG processes
1953 c      print *,"Processor",myrank," fg_rank",fg_rank
1954       if (fg_rank.eq.0) then
1955
1956       if (nodes.eq.1) then
1957         npos=3
1958       else
1959         npos = dlog10(dfloat(nodes-1))+1
1960       endif
1961       if (npos.lt.3) npos=3
1962       write (liczba,'(i1)') npos
1963       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1964      &  //')'
1965       write (liczba,form) me
1966       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1967      &  liczba(:ilen(liczba))
1968       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1969      &  //'.int'
1970       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1971      &  //'.pdb'
1972       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1973      &  liczba(:ilen(liczba))//'.mol2'
1974       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1975      &  liczba(:ilen(liczba))//'.stat'
1976       if (lentmp.gt.0)
1977      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1978      &      //liczba(:ilen(liczba))//'.stat')
1979       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1980      &  //'.rst'
1981       if(usampl) then
1982           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1983      & liczba(:ilen(liczba))//'.const'
1984       endif 
1985
1986       endif
1987 #else
1988       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1989       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1990       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1991       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1992       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1993       if (lentmp.gt.0)
1994      &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
1995      &    //'.stat')
1996       rest2name=prefix(:ilen(prefix))//'.rst'
1997       if(usampl) then 
1998          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1999       endif 
2000 #endif
2001 #if defined(AIX) || defined(PGI)
2002       if (me.eq.king .or. .not. out1file) 
2003      &   open(iout,file=outname,status='unknown')
2004 #ifdef DEBUG
2005       if (fg_rank.gt.0) then
2006         write (liczba,'(i3.3)') myrank/nfgtasks
2007         write (ll,'(bz,i3.3)') fg_rank
2008         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2009      &   status='unknown')
2010       endif
2011 #endif
2012       if(me.eq.king) then
2013        open(igeom,file=intname,status='unknown',position='append')
2014        open(ipdb,file=pdbname,status='unknown')
2015        open(imol2,file=mol2name,status='unknown')
2016        open(istat,file=statname,status='unknown',position='append')
2017       else
2018 c1out       open(iout,file=outname,status='unknown')
2019       endif
2020 #else
2021       if (me.eq.king .or. .not.out1file)
2022      &    open(iout,file=outname,status='unknown')
2023 #ifdef DEBUG
2024       if (fg_rank.gt.0) then
2025         write (liczba,'(i3.3)') myrank/nfgtasks
2026         write (ll,'(bz,i3.3)') fg_rank
2027         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2028      &   status='unknown')
2029       endif
2030 #endif
2031       if(me.eq.king) then
2032        open(igeom,file=intname,status='unknown',access='append')
2033        open(ipdb,file=pdbname,status='unknown')
2034        open(imol2,file=mol2name,status='unknown')
2035        open(istat,file=statname,status='unknown',access='append')
2036       else
2037 c1out       open(iout,file=outname,status='unknown')
2038       endif
2039 #endif
2040       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2041       csa_seed=prefix(:lenpre)//'.CSA.seed'
2042       csa_history=prefix(:lenpre)//'.CSA.history'
2043       csa_bank=prefix(:lenpre)//'.CSA.bank'
2044       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2045       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2046       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2047 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2048       csa_int=prefix(:lenpre)//'.int'
2049       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2050       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2051       csa_in=prefix(:lenpre)//'.CSA.in'
2052 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2053 C Write file names
2054       if (me.eq.king)then
2055       write (iout,'(80(1h-))')
2056       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2057       write (iout,'(80(1h-))')
2058       write (iout,*) "Input file                      : ",
2059      &  pref_orig(:ilen(pref_orig))//'.inp'
2060       write (iout,*) "Output file                     : ",
2061      &  outname(:ilen(outname))
2062       write (iout,*)
2063       write (iout,*) "Sidechain potential file        : ",
2064      &  sidename(:ilen(sidename))
2065 #ifndef OLDSCP
2066       write (iout,*) "SCp potential file              : ",
2067      &  scpname(:ilen(scpname))
2068 #endif
2069       write (iout,*) "Electrostatic potential file    : ",
2070      &  elename(:ilen(elename))
2071       write (iout,*) "Cumulant coefficient file       : ",
2072      &  fouriername(:ilen(fouriername))
2073       write (iout,*) "Torsional parameter file        : ",
2074      &  torname(:ilen(torname))
2075       write (iout,*) "Double torsional parameter file : ",
2076      &  tordname(:ilen(tordname))
2077       write (iout,*) "SCCOR parameter file : ",
2078      &  sccorname(:ilen(sccorname))
2079       write (iout,*) "Bond & inertia constant file    : ",
2080      &  bondname(:ilen(bondname))
2081       write (iout,*) "Bending parameter file          : ",
2082      &  thetname(:ilen(thetname))
2083       write (iout,*) "Rotamer parameter file          : ",
2084      &  rotname(:ilen(rotname))
2085       write (iout,*) "Threading database              : ",
2086      &  patname(:ilen(patname))
2087       if (lentmp.ne.0) 
2088      &write (iout,*)" DIRTMP                          : ",
2089      &  tmpdir(:lentmp)
2090       write (iout,'(80(1h-))')
2091       endif
2092       return
2093       end
2094 c----------------------------------------------------------------------------
2095       subroutine card_concat(card)
2096       implicit real*8 (a-h,o-z)
2097       include 'DIMENSIONS'
2098       include 'COMMON.IOUNITS'
2099       character*(*) card
2100       character*80 karta,ucase
2101       external ilen
2102       read (inp,'(a)') karta
2103       karta=ucase(karta)
2104       card=' '
2105       do while (karta(80:80).eq.'&')
2106         card=card(:ilen(card)+1)//karta(:79)
2107         read (inp,'(a)') karta
2108         karta=ucase(karta)
2109       enddo
2110       card=card(:ilen(card)+1)//karta
2111       return
2112       end
2113 c----------------------------------------------------------------------------------
2114       subroutine readrst
2115       implicit real*8 (a-h,o-z)
2116       include 'DIMENSIONS'
2117       include 'COMMON.CHAIN'
2118       include 'COMMON.IOUNITS'
2119       include 'COMMON.MD'
2120       open(irest2,file=rest2name,status='unknown')
2121       read(irest2,*) totT,EK,potE,totE,t_bath
2122       do i=1,2*nres
2123          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2124       enddo
2125       do i=1,2*nres
2126          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2127       enddo
2128       if(usampl) then
2129              read (irest2,*) iset
2130       endif
2131       close(irest2)
2132       return
2133       end
2134 c---------------------------------------------------------------------------------
2135       subroutine read_fragments
2136       implicit real*8 (a-h,o-z)
2137       include 'DIMENSIONS'
2138 #ifdef MPI
2139       include 'mpif.h'
2140 #endif
2141       include 'COMMON.SETUP'
2142       include 'COMMON.CHAIN'
2143       include 'COMMON.IOUNITS'
2144       include 'COMMON.MD'
2145       include 'COMMON.CONTROL'
2146       read(inp,*) nset,nfrag,npair,nfrag_back
2147       if(me.eq.king.or..not.out1file)
2148      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2149      &  " nfrag_back",nfrag_back
2150       do iset=1,nset
2151          read(inp,*) mset(iset)
2152        do i=1,nfrag
2153          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), 
2154      &     qinfrag(i,iset)
2155          if(me.eq.king.or..not.out1file)
2156      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2157      &     ifrag(2,i,iset), qinfrag(i,iset)
2158        enddo
2159        do i=1,npair
2160         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), 
2161      &    qinpair(i,iset)
2162         if(me.eq.king.or..not.out1file)
2163      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2164      &    ipair(2,i,iset), qinpair(i,iset)
2165        enddo 
2166        do i=1,nfrag_back
2167         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2168      &     wfrag_back(3,i,iset),
2169      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2170         if(me.eq.king.or..not.out1file)
2171      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2172      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2173        enddo 
2174       enddo
2175       return
2176       end
2177 c-------------------------------------------------------------------------------
2178       subroutine read_dist_constr
2179       implicit real*8 (a-h,o-z)
2180       include 'DIMENSIONS'
2181 #ifdef MPI
2182       include 'mpif.h'
2183 #endif
2184       include 'COMMON.SETUP'
2185       include 'COMMON.CONTROL'
2186       include 'COMMON.CHAIN'
2187       include 'COMMON.IOUNITS'
2188       include 'COMMON.SBRIDGE'
2189       integer ifrag_(2,100),ipair_(2,100)
2190       double precision wfrag_(100),wpair_(100)
2191       character*500 controlcard
2192 c      write (iout,*) "Calling read_dist_constr"
2193 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2194 c      call flush(iout)
2195       call card_concat(controlcard)
2196       call readi(controlcard,"NFRAG",nfrag_,0)
2197       call readi(controlcard,"NPAIR",npair_,0)
2198       call readi(controlcard,"NDIST",ndist_,0)
2199       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2200       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2201       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2202       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2203       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2204 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2205 c      write (iout,*) "IFRAG"
2206 c      do i=1,nfrag_
2207 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2208 c      enddo
2209 c      write (iout,*) "IPAIR"
2210 c      do i=1,npair_
2211 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2212 c      enddo
2213       call flush(iout)
2214       do i=1,nfrag_
2215         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2216         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2217      &    ifrag_(2,i)=nstart_sup+nsup-1
2218 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2219         call flush(iout)
2220         if (wfrag_(i).gt.0.0d0) then
2221         do j=ifrag_(1,i),ifrag_(2,i)-1
2222           do k=j+1,ifrag_(2,i)
2223 c            write (iout,*) "j",j," k",k
2224             ddjk=dist(j,k)
2225             if (constr_dist.eq.1) then
2226               nhpb=nhpb+1
2227               ihpb(nhpb)=j
2228               jhpb(nhpb)=k
2229               dhpb(nhpb)=ddjk
2230               forcon(nhpb)=wfrag_(i) 
2231             else if (constr_dist.eq.2) then
2232               if (ddjk.le.dist_cut) then
2233                 nhpb=nhpb+1
2234                 ihpb(nhpb)=j
2235                 jhpb(nhpb)=k
2236                 dhpb(nhpb)=ddjk
2237                 forcon(nhpb)=wfrag_(i) 
2238               endif
2239             else
2240               nhpb=nhpb+1
2241               ihpb(nhpb)=j
2242               jhpb(nhpb)=k
2243               dhpb(nhpb)=ddjk
2244               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2245             endif
2246 #ifdef MPI
2247             if (.not.out1file .or. me.eq.king) 
2248      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2249      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2250 #else
2251             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2252      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2253 #endif
2254           enddo
2255         enddo
2256         endif
2257       enddo
2258       do i=1,npair_
2259         if (wpair_(i).gt.0.0d0) then
2260         ii = ipair_(1,i)
2261         jj = ipair_(2,i)
2262         if (ii.gt.jj) then
2263           itemp=ii
2264           ii=jj
2265           jj=itemp
2266         endif
2267         do j=ifrag_(1,ii),ifrag_(2,ii)
2268           do k=ifrag_(1,jj),ifrag_(2,jj)
2269             nhpb=nhpb+1
2270             ihpb(nhpb)=j
2271             jhpb(nhpb)=k
2272             forcon(nhpb)=wpair_(i)
2273             dhpb(nhpb)=dist(j,k)
2274 #ifdef MPI
2275             if (.not.out1file .or. me.eq.king)
2276      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2277      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2278 #else
2279             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2280      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2281 #endif
2282           enddo
2283         enddo
2284         endif
2285       enddo 
2286       do i=1,ndist_
2287         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2288         if (forcon(nhpb+1).gt.0.0d0) then
2289           nhpb=nhpb+1
2290           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2291 #ifdef MPI
2292           if (.not.out1file .or. me.eq.king)
2293      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2294      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2295 #else
2296           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2297      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2298 #endif
2299         endif
2300       enddo
2301       call flush(iout)
2302       return
2303       end
2304 c-------------------------------------------------------------------------------
2305 #ifdef WINIFL
2306       subroutine flush(iu)
2307       return
2308       end
2309 #endif
2310 #ifdef AIX
2311       subroutine flush(iu)
2312       call flush_(iu)
2313       return
2314       end
2315 #endif
2316 c------------------------------------------------------------------------------
2317       subroutine copy_to_tmp(source)
2318       include "DIMENSIONS"
2319       include "COMMON.IOUNITS"
2320       character*(*) source
2321       character* 256 tmpfile
2322       integer ilen
2323       external ilen
2324       logical ex
2325       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2326       inquire(file=tmpfile,exist=ex)
2327       if (ex) then
2328         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2329      &   " to temporary directory..."
2330         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2331         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2332       endif
2333       return
2334       end
2335 c------------------------------------------------------------------------------
2336       subroutine move_from_tmp(source)
2337       include "DIMENSIONS"
2338       include "COMMON.IOUNITS"
2339       character*(*) source
2340       integer ilen
2341       external ilen
2342       write (*,*) "Moving ",source(:ilen(source)),
2343      & " from temporary directory to working directory"
2344       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2345       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2346       return
2347       end
2348 c------------------------------------------------------------------------------
2349       subroutine random_init(seed)
2350 C
2351 C Initialize random number generator
2352 C
2353       implicit real*8 (a-h,o-z)
2354       include 'DIMENSIONS'
2355 #ifdef MPI
2356       include 'mpif.h'
2357       logical OKRandom, prng_restart
2358       real*8  r1
2359       integer iseed_array(4)
2360 #endif
2361       include 'COMMON.IOUNITS'
2362       include 'COMMON.TIME1'
2363       include 'COMMON.THREAD'
2364       include 'COMMON.SBRIDGE'
2365       include 'COMMON.CONTROL'
2366       include 'COMMON.MCM'
2367       include 'COMMON.MAP'
2368       include 'COMMON.HEADER'
2369       include 'COMMON.CSA'
2370       include 'COMMON.CHAIN'
2371       include 'COMMON.MUCA'
2372       include 'COMMON.MD'
2373       include 'COMMON.FFIELD'
2374       include 'COMMON.SETUP'
2375       iseed=-dint(dabs(seed))
2376       if (iseed.eq.0) then
2377         write (iout,'(/80(1h*)/20x,a/80(1h*))') 
2378      &    'Random seed undefined. The program will stop.'
2379         write (*,'(/80(1h*)/20x,a/80(1h*))') 
2380      &    'Random seed undefined. The program will stop.'
2381 #ifdef MPI
2382         call mpi_finalize(mpi_comm_world,ierr)
2383 #endif
2384         stop 'Bad random seed.'
2385       endif
2386 #ifdef MPI
2387       if (fg_rank.eq.0) then
2388       seed=seed*(me+1)+1
2389 #ifdef AMD64
2390       if(me.eq.king)
2391      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2392       OKRandom = prng_restart(me,iseed)
2393 #else
2394       do i=1,4
2395        tmp=65536.0d0**(4-i)
2396        iseed_array(i) = dint(seed/tmp)
2397        seed=seed-iseed_array(i)*tmp
2398       enddo
2399       if(me.eq.king)
2400      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2401      &                 (iseed_array(i),i=1,4)
2402       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2403      &                 (iseed_array(i),i=1,4)
2404       OKRandom = prng_restart(me,iseed_array)
2405 #endif
2406       if (OKRandom) then
2407 c        r1 = prng_next(me)
2408          r1=ran_number(0.0D0,1.0D0)
2409         if(me.eq.king)
2410      &   write (iout,*) 'ran_num',r1
2411         if (r1.lt.0.0d0) OKRandom=.false.
2412       endif
2413       if (.not.OKRandom) then
2414         write (iout,*) 'PRNG IS NOT WORKING!!!'
2415         print *,'PRNG IS NOT WORKING!!!'
2416         if (me.eq.0) then 
2417          call flush(iout)
2418          call mpi_abort(mpi_comm_world,error_msg,ierr)
2419          stop
2420         else
2421          write (iout,*) 'too many processors for parallel prng'
2422          write (*,*) 'too many processors for parallel prng'
2423          call flush(iout)
2424          stop
2425         endif
2426       endif
2427       endif
2428 #else
2429       call vrndst(iseed)
2430       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2431 #endif
2432       return
2433       end