2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
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
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
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
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'/
85 character*320 controlcard
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
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)
120 safety = 60.0d0*safety
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)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
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 ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
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
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
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"
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
280 if (index(controlcard,'TLIST').gt.0) then
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)
288 if (index(controlcard,'MLIST').gt.0) then
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)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
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
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*400 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 call reada(weightcard,"ATRISS",atriss,0.301D0)
684 call reada(weightcard,"BTRISS",btriss,0.021D0)
685 call reada(weightcard,"CTRISS",ctriss,1.001D0)
686 call reada(weightcard,"DTRISS",dtriss,1.001D0)
687 write (iout,*) "ATRISS=", atriss
688 write (iout,*) "BTRISS=", btriss
689 write (iout,*) "CTRISS=", ctriss
690 write (iout,*) "DTRISS=", dtriss
691 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
693 dyn_ss_mask(i)=.false.
697 dyn_ssbond_ij(i,j)=1.0d300
700 call reada(weightcard,"HT",Ht,0.0D0)
702 ss_depth=ebr/wsc-0.25*eps(1,1)
703 Ht=Ht/wsc-0.25*eps(1,1)
704 akcm=akcm*wstrain/wsc
705 akth=akth*wstrain/wsc
706 akct=akct*wstrain/wsc
707 v1ss=v1ss*wstrain/wsc
708 v2ss=v2ss*wstrain/wsc
709 v3ss=v3ss*wstrain/wsc
711 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
714 if(me.eq.king.or..not.out1file) then
715 write (iout,*) "Parameters of the SS-bond potential:"
716 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
718 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
719 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
720 write (iout,*)" HT",Ht
721 print *,'indpdb=',indpdb,' pdbref=',pdbref
723 if (indpdb.gt.0 .or. pdbref) then
724 read(inp,'(a)') pdbfile
725 if(me.eq.king.or..not.out1file)
726 & write (iout,'(2a)') 'PDB data will be read from file ',
727 & pdbfile(:ilen(pdbfile))
728 open(ipdbin,file=pdbfile,status='old',err=33)
730 33 write (iout,'(a)') 'Error opening PDB file.'
733 c write (iout,*) 'Begin reading pdb data'
736 c write (iout,*) 'Finished reading pdb data'
738 if(me.eq.king.or..not.out1file)
739 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
740 & ' nstart_sup=',nstart_sup
742 itype_pdb(i)=itype(i)
746 nct=nstart_sup+nsup-1
747 call contact(.false.,ncont_ref,icont_ref,co)
750 if(me.eq.king.or..not.out1file)
751 & write(iout,*)'Adding sidechains'
755 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
758 do while (fail.and.nsi.le.maxsi)
759 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
762 if(fail) write(iout,*)'Adding sidechain failed for res ',
763 & i,' after ',nsi,' trials'
768 if (indpdb.eq.0) then
769 C Read sequence if not taken from the pdb file.
771 c print *,'nres=',nres
772 if (iscode.gt.0) then
773 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
775 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
777 C Convert sequence to numeric code
779 itype(i)=rescode(i,sequence(i),iscode)
781 C Assign initial virtual bond lengths
787 vbld(i+nres)=dsc(iabs(itype(i)))
788 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
789 c write (iout,*) "i",i," itype",itype(i),
790 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
794 c print '(20i4)',(itype(i),i=1,nres)
797 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
799 if (itype(i).eq.ntyp1) then
803 else if (iabs(itype(i+1)).ne.20) then
805 else if (iabs(itype(i)).ne.20) then
812 if(me.eq.king.or..not.out1file)then
813 write (iout,*) "ITEL"
815 write (iout,*) i,itype(i),itel(i)
817 print *,'Call Read_Bridge.'
820 C 8/13/98 Set limits to generating the dihedral angles
825 read (inp,*) ndih_constr
826 if (ndih_constr.gt.0) then
828 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
829 if(me.eq.king.or..not.out1file)then
831 & 'There are',ndih_constr,' constraints on phi angles.'
833 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
837 phi0(i)=deg2rad*phi0(i)
838 drange(i)=deg2rad*drange(i)
840 if(me.eq.king.or..not.out1file)
841 & write (iout,*) 'FTORS',ftors
844 phibound(1,ii) = phi0(i)-drange(i)
845 phibound(2,ii) = phi0(i)+drange(i)
852 write (iout,'(a)') 'Boundaries in phi angle sampling:'
854 write (iout,'(a3,i5,2f10.1)')
855 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
861 cd print *,'NNT=',NNT,' NCT=',NCT
862 if (itype(1).eq.ntyp1) nnt=2
863 if (itype(nres).eq.ntyp1) nct=nct-1
865 if(me.eq.king.or..not.out1file)
866 & write (iout,'(a,i3)') 'nsup=',nsup
868 if (nsup.le.(nct-nnt+1)) then
869 do i=0,nct-nnt+1-nsup
870 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
876 & 'Error - sequences to be superposed do not match.'
879 do i=0,nsup-(nct-nnt+1)
880 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
882 nstart_sup=nstart_sup+i
888 & 'Error - sequences to be superposed do not match.'
891 if (nsup.eq.0) nsup=nct-nnt
892 if (nstart_sup.eq.0) nstart_sup=nnt
893 if (nstart_seq.eq.0) nstart_seq=nnt
894 if(me.eq.king.or..not.out1file)
895 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
896 & ' nstart_seq=',nstart_seq
898 c--- Zscore rms -------
899 if (nz_start.eq.0) nz_start=nnt
900 if (nz_end.eq.0 .and. nsup.gt.0) then
902 else if (nz_end.eq.0) then
905 if(me.eq.king.or..not.out1file)then
906 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
907 write (iout,*) 'IZ_SC=',iz_sc
909 c----------------------
912 if (.not.pdbref) then
913 call read_angles(inp,*38)
915 38 write (iout,'(a)') 'Error reading reference structure.'
917 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
918 stop 'Error reading reference structure'
922 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
932 call contact(.true.,ncont_ref,icont_ref,co)
934 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
936 if (constr_dist.gt.0) call read_dist_constr
937 write (iout,*) "After read_dist_constr nhpb",nhpb
939 if(me.eq.king.or..not.out1file)
940 & write (iout,*) 'Contact order:',co
942 if(me.eq.king.or..not.out1file)
943 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
946 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
948 if(me.eq.king.or..not.out1file)
949 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
950 & icont_ref(1,i),' ',
951 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
955 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
956 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
957 & modecalc.ne.10) then
958 C If input structure hasn't been supplied from the PDB file read or generate
960 if (iranconf.eq.0 .and. .not. extconf) then
961 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
962 & write (iout,'(a)') 'Initial geometry will be read in.'
964 read(inp,'(8f10.5)',end=36,err=36)
965 & ((c(l,k),l=1,3),k=1,nres),
966 & ((c(l,k+nres),l=1,3),k=nnt,nct)
967 write (iout,*) "Exit READ_CART"
968 write (iout,'(8f10.5)')
969 & ((c(l,k),l=1,3),k=1,nres),
970 & ((c(l,k+nres),l=1,3),k=nnt,nct)
971 call int_from_cart1(.true.)
972 write (iout,*) "Finish INT_TO_CART"
975 dc(j,i)=c(j,i+1)-c(j,i)
976 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
980 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
982 dc(j,i+nres)=c(j,i+nres)-c(j,i)
983 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
989 call read_angles(inp,*36)
992 36 write (iout,'(a)') 'Error reading angle file.'
994 call mpi_finalize( MPI_COMM_WORLD,IERR )
996 stop 'Error reading angle file.'
998 else if (extconf) then
999 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1000 & write (iout,'(a)') 'Extended chain initial geometry.'
1002 theta(i)=90d0*deg2rad
1005 phi(i)=180d0*deg2rad
1008 alph(i)=110d0*deg2rad
1011 omeg(i)=-120d0*deg2rad
1012 if (itype(i).le.0) omeg(i)=-omeg(i)
1015 if(me.eq.king.or..not.out1file)
1016 & write (iout,'(a)') 'Random-generated initial geometry.'
1020 if (me.eq.king .or. fg_rank.eq.0 .and. (
1021 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1025 call gen_rand_conf(itmp,*30)
1027 30 write (iout,*) 'Failed to generate random conformation',
1028 & ', itrial=',itrial
1029 write (*,*) 'Processor:',me,
1030 & ' Failed to generate random conformation',
1040 write (iout,'(a,i3,a)') 'Processor:',me,
1041 & ' error in generating random conformation.'
1042 write (*,'(a,i3,a)') 'Processor:',me,
1043 & ' error in generating random conformation.'
1046 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1052 call gen_rand_conf(itmp,*30)
1054 30 write (iout,*) 'Failed to generate random conformation',
1055 & ', itrial=',itrial
1056 write (*,*) 'Failed to generate random conformation',
1057 & ', itrial=',itrial
1059 write (iout,'(a,i3,a)') 'Processor:',me,
1060 & ' error in generating random conformation.'
1061 write (*,'(a,i3,a)') 'Processor:',me,
1062 & ' error in generating random conformation.'
1067 elseif (modecalc.eq.4) then
1068 read (inp,'(a)') intinname
1069 open (intin,file=intinname,status='old',err=333)
1070 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1071 & write (iout,'(a)') 'intinname',intinname
1072 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1074 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1076 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1078 stop 'Error opening angle file.'
1082 C Generate distance constraints, if the PDB structure is to be regularized.
1083 if (nthread.gt.0) then
1084 call read_threadbase
1087 if (me.eq.king .or. .not. out1file)
1089 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1090 write (iout,'(/a,i3,a)')
1091 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1092 write (iout,'(20i4)') (iss(i),i=1,ns)
1094 write(iout,*)"Running with dynamic disulfide-bond formation"
1096 write (iout,'(/a/)') 'Pre-formed links are:'
1102 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1103 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1109 if (ns.gt.0.and.dyn_ss) then
1113 forcon(i-nss)=forcon(i)
1120 dyn_ss_mask(iss(i))=.true.
1123 if (i2ndstr.gt.0) call secstrp2dihc
1124 c call geom_to_var(nvar,x)
1125 c call etotal(energia(0))
1126 c call enerprint(energia(0))
1127 c call briefout(0,etot)
1129 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1130 cd write (iout,'(a)') 'Variable list:'
1131 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1133 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1134 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1135 & 'Processor',myrank,': end reading molecular data.'
1139 c--------------------------------------------------------------------------
1140 logical function seq_comp(itypea,itypeb,length)
1142 integer length,itypea(length),itypeb(length)
1145 if (itypea(i).ne.itypeb(i)) then
1153 c-----------------------------------------------------------------------------
1154 subroutine read_bridge
1155 C Read information about disulfide bridges.
1156 implicit real*8 (a-h,o-z)
1157 include 'DIMENSIONS'
1161 include 'COMMON.IOUNITS'
1162 include 'COMMON.GEO'
1163 include 'COMMON.VAR'
1164 include 'COMMON.INTERACT'
1165 include 'COMMON.LOCAL'
1166 include 'COMMON.NAMES'
1167 include 'COMMON.CHAIN'
1168 include 'COMMON.FFIELD'
1169 include 'COMMON.SBRIDGE'
1170 include 'COMMON.HEADER'
1171 include 'COMMON.CONTROL'
1172 include 'COMMON.DBASE'
1173 include 'COMMON.THREAD'
1174 include 'COMMON.TIME1'
1175 include 'COMMON.SETUP'
1176 C Read bridging residues.
1177 read (inp,*) ns,(iss(i),i=1,ns)
1179 if(me.eq.king.or..not.out1file)
1180 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1181 C Check whether the specified bridging residues are cystines.
1183 if (itype(iss(i)).ne.1) then
1184 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1185 & 'Do you REALLY think that the residue ',
1186 & restyp(itype(iss(i))),i,
1187 & ' can form a disulfide bridge?!!!'
1188 write (*,'(2a,i3,a)')
1189 & 'Do you REALLY think that the residue ',
1190 & restyp(itype(iss(i))),i,
1191 & ' can form a disulfide bridge?!!!'
1193 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1198 C Read preformed bridges.
1200 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1202 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1205 C Check if the residues involved in bridges are in the specified list of
1206 C bridging residues.
1209 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1210 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1211 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1212 & ' contains residues present in other pairs.'
1213 write (*,'(a,i3,a)') 'Disulfide pair',i,
1214 & ' contains residues present in other pairs.'
1216 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1222 if (ihpb(i).eq.iss(j)) goto 10
1224 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1227 if (jhpb(i).eq.iss(j)) goto 20
1229 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1235 ihpb(i)=ihpb(i)+nres
1236 jhpb(i)=jhpb(i)+nres
1242 c----------------------------------------------------------------------------
1243 subroutine read_x(kanal,*)
1244 implicit real*8 (a-h,o-z)
1245 include 'DIMENSIONS'
1246 include 'COMMON.GEO'
1247 include 'COMMON.VAR'
1248 include 'COMMON.CHAIN'
1249 include 'COMMON.IOUNITS'
1250 include 'COMMON.CONTROL'
1251 include 'COMMON.LOCAL'
1252 include 'COMMON.INTERACT'
1253 c Read coordinates from input
1255 read(kanal,'(8f10.5)',end=10,err=10)
1256 & ((c(l,k),l=1,3),k=1,nres),
1257 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1260 c(j,2*nres)=c(j,nres)
1262 call int_from_cart1(.false.)
1265 dc(j,i)=c(j,i+1)-c(j,i)
1266 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1270 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1272 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1273 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1281 c----------------------------------------------------------------------------
1282 subroutine read_threadbase
1283 implicit real*8 (a-h,o-z)
1284 include 'DIMENSIONS'
1285 include 'COMMON.IOUNITS'
1286 include 'COMMON.GEO'
1287 include 'COMMON.VAR'
1288 include 'COMMON.INTERACT'
1289 include 'COMMON.LOCAL'
1290 include 'COMMON.NAMES'
1291 include 'COMMON.CHAIN'
1292 include 'COMMON.FFIELD'
1293 include 'COMMON.SBRIDGE'
1294 include 'COMMON.HEADER'
1295 include 'COMMON.CONTROL'
1296 include 'COMMON.DBASE'
1297 include 'COMMON.THREAD'
1298 include 'COMMON.TIME1'
1299 C Read pattern database for threading.
1300 read (icbase,*) nseq
1302 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1303 & nres_base(2,i),nres_base(3,i)
1304 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1306 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1307 c & nres_base(2,i),nres_base(3,i)
1308 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1312 if (weidis.eq.0.0D0) weidis=0.1D0
1321 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1322 write (iout,'(a,i5)') 'nexcl: ',nexcl
1323 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1326 c------------------------------------------------------------------------------
1327 subroutine setup_var
1328 implicit real*8 (a-h,o-z)
1329 include 'DIMENSIONS'
1330 include 'COMMON.IOUNITS'
1331 include 'COMMON.GEO'
1332 include 'COMMON.VAR'
1333 include 'COMMON.INTERACT'
1334 include 'COMMON.LOCAL'
1335 include 'COMMON.NAMES'
1336 include 'COMMON.CHAIN'
1337 include 'COMMON.FFIELD'
1338 include 'COMMON.SBRIDGE'
1339 include 'COMMON.HEADER'
1340 include 'COMMON.CONTROL'
1341 include 'COMMON.DBASE'
1342 include 'COMMON.THREAD'
1343 include 'COMMON.TIME1'
1344 C Set up variable list.
1350 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1352 ialph(i,1)=nvar+nside
1356 if (indphi.gt.0) then
1358 else if (indback.gt.0) then
1363 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1366 c----------------------------------------------------------------------------
1367 subroutine gen_dist_constr
1368 C Generate CA distance constraints.
1369 implicit real*8 (a-h,o-z)
1370 include 'DIMENSIONS'
1371 include 'COMMON.IOUNITS'
1372 include 'COMMON.GEO'
1373 include 'COMMON.VAR'
1374 include 'COMMON.INTERACT'
1375 include 'COMMON.LOCAL'
1376 include 'COMMON.NAMES'
1377 include 'COMMON.CHAIN'
1378 include 'COMMON.FFIELD'
1379 include 'COMMON.SBRIDGE'
1380 include 'COMMON.HEADER'
1381 include 'COMMON.CONTROL'
1382 include 'COMMON.DBASE'
1383 include 'COMMON.THREAD'
1384 include 'COMMON.TIME1'
1385 dimension itype_pdb(maxres)
1386 common /pizda/ itype_pdb
1388 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1389 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1390 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1392 do i=nstart_sup,nstart_sup+nsup-1
1393 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1394 cd & ' seq_pdb', restyp(itype_pdb(i))
1395 do j=i+2,nstart_sup+nsup-1
1397 ihpb(nhpb)=i+nstart_seq-nstart_sup
1398 jhpb(nhpb)=j+nstart_seq-nstart_sup
1400 dhpb(nhpb)=dist(i,j)
1403 cd write (iout,'(a)') 'Distance constraints:'
1408 cd if (ii.gt.nres) then
1413 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1414 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1415 cd & dhpb(i),forcon(i)
1419 c----------------------------------------------------------------------------
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1423 include 'COMMON.MAP'
1424 include 'COMMON.IOUNITS'
1425 character*3 angid(4) /'THE','PHI','ALP','OME'/
1426 character*80 mapcard,ucase
1428 read (inp,'(a)') mapcard
1429 mapcard=ucase(mapcard)
1430 if (index(mapcard,'PHI').gt.0) then
1432 else if (index(mapcard,'THE').gt.0) then
1434 else if (index(mapcard,'ALP').gt.0) then
1436 else if (index(mapcard,'OME').gt.0) then
1439 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1440 stop 'Error - illegal variable spec in MAP card.'
1442 call readi (mapcard,'RES1',res1(imap),0)
1443 call readi (mapcard,'RES2',res2(imap),0)
1444 if (res1(imap).eq.0) then
1445 res1(imap)=res2(imap)
1446 else if (res2(imap).eq.0) then
1447 res2(imap)=res1(imap)
1449 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1451 & 'Error - illegal definition of variable group in MAP.'
1452 stop 'Error - illegal definition of variable group in MAP.'
1454 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1455 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1456 call readi(mapcard,'NSTEP',nstep(imap),0)
1457 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1459 & 'Illegal boundary and/or step size specification in MAP.'
1460 stop 'Illegal boundary and/or step size specification in MAP.'
1465 c----------------------------------------------------------------------------
1467 implicit real*8 (a-h,o-z)
1468 include 'DIMENSIONS'
1469 include 'COMMON.IOUNITS'
1470 include 'COMMON.GEO'
1471 include 'COMMON.CSA'
1472 include 'COMMON.BANK'
1473 include 'COMMON.CONTROL'
1475 character*620 mcmcard
1476 call card_concat(mcmcard)
1478 call readi(mcmcard,'NCONF',nconf,50)
1479 call readi(mcmcard,'NADD',nadd,0)
1480 call readi(mcmcard,'JSTART',jstart,1)
1481 call readi(mcmcard,'JEND',jend,1)
1482 call readi(mcmcard,'NSTMAX',nstmax,500000)
1483 call readi(mcmcard,'N0',n0,1)
1484 call readi(mcmcard,'N1',n1,6)
1485 call readi(mcmcard,'N2',n2,4)
1486 call readi(mcmcard,'N3',n3,0)
1487 call readi(mcmcard,'N4',n4,0)
1488 call readi(mcmcard,'N5',n5,0)
1489 call readi(mcmcard,'N6',n6,10)
1490 call readi(mcmcard,'N7',n7,0)
1491 call readi(mcmcard,'N8',n8,0)
1492 call readi(mcmcard,'N9',n9,0)
1493 call readi(mcmcard,'N14',n14,0)
1494 call readi(mcmcard,'N15',n15,0)
1495 call readi(mcmcard,'N16',n16,0)
1496 call readi(mcmcard,'N17',n17,0)
1497 call readi(mcmcard,'N18',n18,0)
1499 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1501 call readi(mcmcard,'NDIFF',ndiff,2)
1502 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1503 call readi(mcmcard,'IS1',is1,1)
1504 call readi(mcmcard,'IS2',is2,8)
1505 call readi(mcmcard,'NRAN0',nran0,4)
1506 call readi(mcmcard,'NRAN1',nran1,2)
1507 call readi(mcmcard,'IRR',irr,1)
1508 call readi(mcmcard,'NSEED',nseed,20)
1509 call readi(mcmcard,'NTOTAL',ntotal,10000)
1510 call reada(mcmcard,'CUT1',cut1,2.0d0)
1511 call reada(mcmcard,'CUT2',cut2,5.0d0)
1512 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1513 call readi(mcmcard,'ICMAX',icmax,3)
1514 call readi(mcmcard,'IRESTART',irestart,0)
1515 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1518 call reada(mcmcard,'DELE',dele,20.0d0)
1519 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1520 call readi(mcmcard,'IREF',iref,0)
1521 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1522 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1523 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1524 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1525 write (iout,*) "NCONF_IN",nconf_in
1528 c----------------------------------------------------------------------------
1529 cfmc subroutine mcmfread
1530 cfmc implicit real*8 (a-h,o-z)
1531 cfmc include 'DIMENSIONS'
1532 cfmc include 'COMMON.MCMF'
1533 cfmc include 'COMMON.IOUNITS'
1534 cfmc include 'COMMON.GEO'
1535 cfmc character*80 ucase
1536 cfmc character*620 mcmcard
1537 cfmc call card_concat(mcmcard)
1539 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1540 cfmc write(iout,*)'MAXRANT=',maxrant
1541 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1542 cfmc write(iout,*)'MAXFAM=',maxfam
1543 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1544 cfmc write(iout,*)'NNET1=',nnet1
1545 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1546 cfmc write(iout,*)'NNET2=',nnet2
1547 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1548 cfmc write(iout,*)'NNET3=',nnet3
1549 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1550 cfmc write(iout,*)'ILASTT=',ilastt
1551 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1552 cfmc write(iout,*)'MAXSTR=',maxstr
1553 cfmc maxstr_f=maxstr/maxfam
1554 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1555 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1556 cfmc write(iout,*)'NMCMF=',nmcmf
1557 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1558 cfmc write(iout,*)'IFOCUS=',ifocus
1559 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1560 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1561 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1562 cfmc write(iout,*)'INTPRT=',intprt
1563 cfmc call readi(mcmcard,'IPRT',iprt,100)
1564 cfmc write(iout,*)'IPRT=',iprt
1565 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1566 cfmc write(iout,*)'IMAXTR=',imaxtr
1567 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1568 cfmc write(iout,*)'MAXEVEN=',maxeven
1569 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1570 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1571 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1572 cfmc write(iout,*)'INIMIN=',inimin
1573 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1574 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1575 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1576 cfmc write(iout,*)'NTHREAD=',nthread
1577 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1578 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1579 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1580 cfmc write(iout,*)'MAXPERT=',maxpert
1581 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1582 cfmc write(iout,*)'IRMSD=',irmsd
1583 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1584 cfmc write(iout,*)'DENEMIN=',denemin
1585 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1586 cfmc write(iout,*)'RCUT1S=',rcut1s
1587 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1588 cfmc write(iout,*)'RCUT1E=',rcut1e
1589 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1590 cfmc write(iout,*)'RCUT2S=',rcut2s
1591 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1592 cfmc write(iout,*)'RCUT2E=',rcut2e
1593 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1594 cfmc write(iout,*)'DPERT1=',d_pert1
1595 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1596 cfmc write(iout,*)'DPERT1A=',d_pert1a
1597 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1598 cfmc write(iout,*)'DPERT2=',d_pert2
1599 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1600 cfmc write(iout,*)'DPERT2A=',d_pert2a
1601 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1602 cfmc write(iout,*)'DPERT2B=',d_pert2b
1603 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1604 cfmc write(iout,*)'DPERT2C=',d_pert2c
1605 cfmc d_pert1=deg2rad*d_pert1
1606 cfmc d_pert1a=deg2rad*d_pert1a
1607 cfmc d_pert2=deg2rad*d_pert2
1608 cfmc d_pert2a=deg2rad*d_pert2a
1609 cfmc d_pert2b=deg2rad*d_pert2b
1610 cfmc d_pert2c=deg2rad*d_pert2c
1611 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1612 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1613 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1614 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1615 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1616 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1617 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1618 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1619 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1620 cfmc write(iout,*)'RCUTINI=',rcutini
1621 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1622 cfmc write(iout,*)'GRAT=',grat
1623 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1624 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1628 c----------------------------------------------------------------------------
1630 implicit real*8 (a-h,o-z)
1631 include 'DIMENSIONS'
1632 include 'COMMON.MCM'
1633 include 'COMMON.MCE'
1634 include 'COMMON.IOUNITS'
1636 character*320 mcmcard
1637 call card_concat(mcmcard)
1638 call readi(mcmcard,'MAXACC',maxacc,100)
1639 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1640 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1641 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1642 call readi(mcmcard,'MAXREPM',maxrepm,200)
1643 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1644 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1645 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1646 call reada(mcmcard,'E_UP',e_up,5.0D0)
1647 call reada(mcmcard,'DELTE',delte,0.1D0)
1648 call readi(mcmcard,'NSWEEP',nsweep,5)
1649 call readi(mcmcard,'NSTEPH',nsteph,0)
1650 call readi(mcmcard,'NSTEPC',nstepc,0)
1651 call reada(mcmcard,'TMIN',tmin,298.0D0)
1652 call reada(mcmcard,'TMAX',tmax,298.0D0)
1653 call readi(mcmcard,'NWINDOW',nwindow,0)
1654 call readi(mcmcard,'PRINT_MC',print_mc,0)
1655 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1656 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1657 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1658 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1659 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1660 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1661 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1662 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1663 if (nwindow.gt.0) then
1664 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1666 winlen(i)=winend(i)-winstart(i)+1
1669 if (tmax.lt.tmin) tmax=tmin
1670 if (tmax.eq.tmin) then
1674 if (nstepc.gt.0 .and. nsteph.gt.0) then
1675 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1676 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1678 C Probabilities of different move types
1679 sumpro_type(0)=0.0D0
1680 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1681 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1682 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1683 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1684 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1685 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1686 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1688 print *,'i',i,' sumprotype',sumpro_type(i)
1689 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1690 print *,'i',i,' sumprotype',sumpro_type(i)
1694 c----------------------------------------------------------------------------
1695 subroutine read_minim
1696 implicit real*8 (a-h,o-z)
1697 include 'DIMENSIONS'
1698 include 'COMMON.MINIM'
1699 include 'COMMON.IOUNITS'
1701 character*320 minimcard
1702 call card_concat(minimcard)
1703 call readi(minimcard,'MAXMIN',maxmin,2000)
1704 call readi(minimcard,'MAXFUN',maxfun,5000)
1705 call readi(minimcard,'MINMIN',minmin,maxmin)
1706 call readi(minimcard,'MINFUN',minfun,maxmin)
1707 call reada(minimcard,'TOLF',tolf,1.0D-2)
1708 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1709 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1710 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1711 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1712 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1713 & 'Options in energy minimization:'
1714 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1715 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1716 & 'MinMin:',MinMin,' MinFun:',MinFun,
1717 & ' TolF:',TolF,' RTolF:',RTolF
1720 c----------------------------------------------------------------------------
1721 subroutine read_angles(kanal,*)
1722 implicit real*8 (a-h,o-z)
1723 include 'DIMENSIONS'
1724 include 'COMMON.GEO'
1725 include 'COMMON.VAR'
1726 include 'COMMON.CHAIN'
1727 include 'COMMON.IOUNITS'
1728 include 'COMMON.CONTROL'
1729 c Read angles from input
1731 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1732 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1733 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1734 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1737 c 9/7/01 avoid 180 deg valence angle
1738 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1740 theta(i)=deg2rad*theta(i)
1741 phi(i)=deg2rad*phi(i)
1742 alph(i)=deg2rad*alph(i)
1743 omeg(i)=deg2rad*omeg(i)
1748 c----------------------------------------------------------------------------
1749 subroutine reada(rekord,lancuch,wartosc,default)
1751 character*(*) rekord,lancuch
1752 double precision wartosc,default
1755 iread=index(rekord,lancuch)
1756 if (iread.eq.0) then
1760 iread=iread+ilen(lancuch)+1
1761 read (rekord(iread:),*,err=10,end=10) wartosc
1766 c----------------------------------------------------------------------------
1767 subroutine readi(rekord,lancuch,wartosc,default)
1769 character*(*) rekord,lancuch
1770 integer wartosc,default
1773 iread=index(rekord,lancuch)
1774 if (iread.eq.0) then
1778 iread=iread+ilen(lancuch)+1
1779 read (rekord(iread:),*,err=10,end=10) wartosc
1784 c----------------------------------------------------------------------------
1785 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1788 integer tablica(dim),default
1789 character*(*) rekord,lancuch
1796 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1797 if (iread.eq.0) return
1798 iread=iread+ilen(lancuch)+1
1799 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1802 c----------------------------------------------------------------------------
1803 subroutine multreada(rekord,lancuch,tablica,dim,default)
1806 double precision tablica(dim),default
1807 character*(*) rekord,lancuch
1814 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1815 if (iread.eq.0) return
1816 iread=iread+ilen(lancuch)+1
1817 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1820 c----------------------------------------------------------------------------
1821 subroutine openunits
1822 implicit real*8 (a-h,o-z)
1823 include 'DIMENSIONS'
1826 character*16 form,nodename
1829 include 'COMMON.SETUP'
1830 include 'COMMON.IOUNITS'
1832 include 'COMMON.CONTROL'
1833 integer lenpre,lenpot,ilen,lentmp
1835 character*3 out1file_text,ucase
1838 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1839 call getenv_loc("PREFIX",prefix)
1841 call getenv_loc("POT",pot)
1842 call getenv_loc("DIRTMP",tmpdir)
1843 call getenv_loc("CURDIR",curdir)
1844 call getenv_loc("OUT1FILE",out1file_text)
1845 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1846 out1file_text=ucase(out1file_text)
1847 if (out1file_text(1:1).eq."Y") then
1850 out1file=fg_rank.gt.0
1855 if (lentmp.gt.0) then
1856 write (*,'(80(1h!))')
1857 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1858 write (*,'(80(1h!))')
1859 write (*,*)"All output files will be on node /tmp directory."
1861 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1862 if (me.eq.king) then
1863 write (*,*) "The master node is ",nodename
1864 else if (fg_rank.eq.0) then
1865 write (*,*) "I am the CG slave node ",nodename
1867 write (*,*) "I am the FG slave node ",nodename
1870 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1871 lenpre = lentmp+lenpre+1
1873 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1874 C Get the names and open the input files
1875 #if defined(WINIFL) || defined(WINPGI)
1876 open(1,file=pref_orig(:ilen(pref_orig))//
1877 & '.inp',status='old',readonly,shared)
1878 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1879 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1880 C Get parameter filenames and open the parameter files.
1881 call getenv_loc('BONDPAR',bondname)
1882 open (ibond,file=bondname,status='old',readonly,shared)
1883 call getenv_loc('THETPAR',thetname)
1884 open (ithep,file=thetname,status='old',readonly,shared)
1885 call getenv_loc('ROTPAR',rotname)
1886 open (irotam,file=rotname,status='old',readonly,shared)
1887 call getenv_loc('TORPAR',torname)
1888 open (itorp,file=torname,status='old',readonly,shared)
1889 call getenv_loc('TORDPAR',tordname)
1890 open (itordp,file=tordname,status='old',readonly,shared)
1891 call getenv_loc('FOURIER',fouriername)
1892 open (ifourier,file=fouriername,status='old',readonly,shared)
1893 call getenv_loc('ELEPAR',elename)
1894 open (ielep,file=elename,status='old',readonly,shared)
1895 call getenv_loc('SIDEPAR',sidename)
1896 open (isidep,file=sidename,status='old',readonly,shared)
1897 #elif (defined CRAY) || (defined AIX)
1898 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1900 c print *,"Processor",myrank," opened file 1"
1901 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1902 c print *,"Processor",myrank," opened file 9"
1903 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1904 C Get parameter filenames and open the parameter files.
1905 call getenv_loc('BONDPAR',bondname)
1906 open (ibond,file=bondname,status='old',action='read')
1907 c print *,"Processor",myrank," opened file IBOND"
1908 call getenv_loc('THETPAR',thetname)
1909 open (ithep,file=thetname,status='old',action='read')
1910 c print *,"Processor",myrank," opened file ITHEP"
1911 call getenv_loc('ROTPAR',rotname)
1912 open (irotam,file=rotname,status='old',action='read')
1913 c print *,"Processor",myrank," opened file IROTAM"
1914 call getenv_loc('TORPAR',torname)
1915 open (itorp,file=torname,status='old',action='read')
1916 c print *,"Processor",myrank," opened file ITORP"
1917 call getenv_loc('TORDPAR',tordname)
1918 open (itordp,file=tordname,status='old',action='read')
1919 c print *,"Processor",myrank," opened file ITORDP"
1920 call getenv_loc('SCCORPAR',sccorname)
1921 open (isccor,file=sccorname,status='old',action='read')
1922 c print *,"Processor",myrank," opened file ISCCOR"
1923 call getenv_loc('FOURIER',fouriername)
1924 open (ifourier,file=fouriername,status='old',action='read')
1925 c print *,"Processor",myrank," opened file IFOURIER"
1926 call getenv_loc('ELEPAR',elename)
1927 open (ielep,file=elename,status='old',action='read')
1928 c print *,"Processor",myrank," opened file IELEP"
1929 call getenv_loc('SIDEPAR',sidename)
1930 open (isidep,file=sidename,status='old',action='read')
1931 c print *,"Processor",myrank," opened file ISIDEP"
1932 c print *,"Processor",myrank," opened parameter files"
1934 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1935 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1936 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1937 C Get parameter filenames and open the parameter files.
1938 call getenv_loc('BONDPAR',bondname)
1939 open (ibond,file=bondname,status='old')
1940 call getenv_loc('THETPAR',thetname)
1941 open (ithep,file=thetname,status='old')
1942 call getenv_loc('ROTPAR',rotname)
1943 open (irotam,file=rotname,status='old')
1944 call getenv_loc('TORPAR',torname)
1945 open (itorp,file=torname,status='old')
1946 call getenv_loc('TORDPAR',tordname)
1947 open (itordp,file=tordname,status='old')
1948 call getenv_loc('SCCORPAR',sccorname)
1949 open (isccor,file=sccorname,status='old')
1950 call getenv_loc('FOURIER',fouriername)
1951 open (ifourier,file=fouriername,status='old')
1952 call getenv_loc('ELEPAR',elename)
1953 open (ielep,file=elename,status='old')
1954 call getenv_loc('SIDEPAR',sidename)
1955 open (isidep,file=sidename,status='old')
1957 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1959 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1960 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1961 C Get parameter filenames and open the parameter files.
1962 call getenv_loc('BONDPAR',bondname)
1963 open (ibond,file=bondname,status='old',readonly)
1964 call getenv_loc('THETPAR',thetname)
1965 open (ithep,file=thetname,status='old',readonly)
1966 call getenv_loc('ROTPAR',rotname)
1967 open (irotam,file=rotname,status='old',readonly)
1968 call getenv_loc('TORPAR',torname)
1969 open (itorp,file=torname,status='old',readonly)
1970 call getenv_loc('TORDPAR',tordname)
1971 open (itordp,file=tordname,status='old',readonly)
1972 call getenv_loc('SCCORPAR',sccorname)
1973 open (isccor,file=sccorname,status='old',readonly)
1975 call getenv_loc('THETPARPDB',thetname_pdb)
1976 print *,"thetname_pdb ",thetname_pdb
1977 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1978 print *,ithep_pdb," opened"
1980 call getenv_loc('FOURIER',fouriername)
1981 open (ifourier,file=fouriername,status='old',readonly)
1982 call getenv_loc('ELEPAR',elename)
1983 open (ielep,file=elename,status='old',readonly)
1984 call getenv_loc('SIDEPAR',sidename)
1985 open (isidep,file=sidename,status='old',readonly)
1987 call getenv_loc('ROTPARPDB',rotname_pdb)
1988 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1993 C 8/9/01 In the newest version SCp interaction constants are read from a file
1994 C Use -DOLDSCP to use hard-coded constants instead.
1996 call getenv_loc('SCPPAR',scpname)
1997 #if defined(WINIFL) || defined(WINPGI)
1998 open (iscpp,file=scpname,status='old',readonly,shared)
1999 #elif (defined CRAY) || (defined AIX)
2000 open (iscpp,file=scpname,status='old',action='read')
2002 open (iscpp,file=scpname,status='old')
2004 open (iscpp,file=scpname,status='old',readonly)
2007 call getenv_loc('PATTERN',patname)
2008 #if defined(WINIFL) || defined(WINPGI)
2009 open (icbase,file=patname,status='old',readonly,shared)
2010 #elif (defined CRAY) || (defined AIX)
2011 open (icbase,file=patname,status='old',action='read')
2013 open (icbase,file=patname,status='old')
2015 open (icbase,file=patname,status='old',readonly)
2018 C Open output file only for CG processes
2019 c print *,"Processor",myrank," fg_rank",fg_rank
2020 if (fg_rank.eq.0) then
2022 if (nodes.eq.1) then
2025 npos = dlog10(dfloat(nodes-1))+1
2027 if (npos.lt.3) npos=3
2028 write (liczba,'(i1)') npos
2029 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2031 write (liczba,form) me
2032 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2033 & liczba(:ilen(liczba))
2034 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2036 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2038 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2039 & liczba(:ilen(liczba))//'.mol2'
2040 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2041 & liczba(:ilen(liczba))//'.stat'
2043 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2044 & //liczba(:ilen(liczba))//'.stat')
2045 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2048 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2049 & liczba(:ilen(liczba))//'.const'
2054 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2055 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2056 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2057 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2058 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2060 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2062 rest2name=prefix(:ilen(prefix))//'.rst'
2064 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2067 #if defined(AIX) || defined(PGI)
2068 if (me.eq.king .or. .not. out1file)
2069 & open(iout,file=outname,status='unknown')
2071 if (fg_rank.gt.0) then
2072 write (liczba,'(i3.3)') myrank/nfgtasks
2073 write (ll,'(bz,i3.3)') fg_rank
2074 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2079 open(igeom,file=intname,status='unknown',position='append')
2080 open(ipdb,file=pdbname,status='unknown')
2081 open(imol2,file=mol2name,status='unknown')
2082 open(istat,file=statname,status='unknown',position='append')
2084 c1out open(iout,file=outname,status='unknown')
2087 if (me.eq.king .or. .not.out1file)
2088 & open(iout,file=outname,status='unknown')
2090 if (fg_rank.gt.0) then
2091 write (liczba,'(i3.3)') myrank/nfgtasks
2092 write (ll,'(bz,i3.3)') fg_rank
2093 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2098 open(igeom,file=intname,status='unknown',access='append')
2099 open(ipdb,file=pdbname,status='unknown')
2100 open(imol2,file=mol2name,status='unknown')
2101 open(istat,file=statname,status='unknown',access='append')
2103 c1out open(iout,file=outname,status='unknown')
2106 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2107 csa_seed=prefix(:lenpre)//'.CSA.seed'
2108 csa_history=prefix(:lenpre)//'.CSA.history'
2109 csa_bank=prefix(:lenpre)//'.CSA.bank'
2110 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2111 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2112 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2113 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2114 csa_int=prefix(:lenpre)//'.int'
2115 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2116 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2117 csa_in=prefix(:lenpre)//'.CSA.in'
2118 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2121 write (iout,'(80(1h-))')
2122 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2123 write (iout,'(80(1h-))')
2124 write (iout,*) "Input file : ",
2125 & pref_orig(:ilen(pref_orig))//'.inp'
2126 write (iout,*) "Output file : ",
2127 & outname(:ilen(outname))
2129 write (iout,*) "Sidechain potential file : ",
2130 & sidename(:ilen(sidename))
2132 write (iout,*) "SCp potential file : ",
2133 & scpname(:ilen(scpname))
2135 write (iout,*) "Electrostatic potential file : ",
2136 & elename(:ilen(elename))
2137 write (iout,*) "Cumulant coefficient file : ",
2138 & fouriername(:ilen(fouriername))
2139 write (iout,*) "Torsional parameter file : ",
2140 & torname(:ilen(torname))
2141 write (iout,*) "Double torsional parameter file : ",
2142 & tordname(:ilen(tordname))
2143 write (iout,*) "SCCOR parameter file : ",
2144 & sccorname(:ilen(sccorname))
2145 write (iout,*) "Bond & inertia constant file : ",
2146 & bondname(:ilen(bondname))
2147 write (iout,*) "Bending parameter file : ",
2148 & thetname(:ilen(thetname))
2149 write (iout,*) "Rotamer parameter file : ",
2150 & rotname(:ilen(rotname))
2151 write (iout,*) "Thetpdb parameter file : ",
2152 & thetname_pdb(:ilen(thetname_pdb))
2153 write (iout,*) "Threading database : ",
2154 & patname(:ilen(patname))
2156 &write (iout,*)" DIRTMP : ",
2158 write (iout,'(80(1h-))')
2162 c----------------------------------------------------------------------------
2163 subroutine card_concat(card)
2164 implicit real*8 (a-h,o-z)
2165 include 'DIMENSIONS'
2166 include 'COMMON.IOUNITS'
2168 character*80 karta,ucase
2170 read (inp,'(a)') karta
2173 do while (karta(80:80).eq.'&')
2174 card=card(:ilen(card)+1)//karta(:79)
2175 read (inp,'(a)') karta
2178 card=card(:ilen(card)+1)//karta
2181 c----------------------------------------------------------------------------------
2183 implicit real*8 (a-h,o-z)
2184 include 'DIMENSIONS'
2185 include 'COMMON.CHAIN'
2186 include 'COMMON.IOUNITS'
2188 open(irest2,file=rest2name,status='unknown')
2189 read(irest2,*) totT,EK,potE,totE,t_bath
2191 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2194 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2197 read (irest2,*) iset
2202 c---------------------------------------------------------------------------------
2203 subroutine read_fragments
2204 implicit real*8 (a-h,o-z)
2205 include 'DIMENSIONS'
2209 include 'COMMON.SETUP'
2210 include 'COMMON.CHAIN'
2211 include 'COMMON.IOUNITS'
2213 include 'COMMON.CONTROL'
2214 read(inp,*) nset,nfrag,npair,nfrag_back
2215 if(me.eq.king.or..not.out1file)
2216 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2217 & " nfrag_back",nfrag_back
2219 read(inp,*) mset(iset)
2221 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2223 if(me.eq.king.or..not.out1file)
2224 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2225 & ifrag(2,i,iset), qinfrag(i,iset)
2228 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2230 if(me.eq.king.or..not.out1file)
2231 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2232 & ipair(2,i,iset), qinpair(i,iset)
2235 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2236 & wfrag_back(3,i,iset),
2237 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2238 if(me.eq.king.or..not.out1file)
2239 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2240 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2245 c-------------------------------------------------------------------------------
2246 subroutine read_dist_constr
2247 implicit real*8 (a-h,o-z)
2248 include 'DIMENSIONS'
2252 include 'COMMON.SETUP'
2253 include 'COMMON.CONTROL'
2254 include 'COMMON.CHAIN'
2255 include 'COMMON.IOUNITS'
2256 include 'COMMON.SBRIDGE'
2257 integer ifrag_(2,100),ipair_(2,100)
2258 double precision wfrag_(100),wpair_(100)
2259 character*500 controlcard
2260 c write (iout,*) "Calling read_dist_constr"
2261 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2263 call card_concat(controlcard)
2264 call readi(controlcard,"NFRAG",nfrag_,0)
2265 call readi(controlcard,"NPAIR",npair_,0)
2266 call readi(controlcard,"NDIST",ndist_,0)
2267 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2268 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2269 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2270 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2271 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2272 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2273 c write (iout,*) "IFRAG"
2275 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2277 c write (iout,*) "IPAIR"
2279 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2283 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2284 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2285 & ifrag_(2,i)=nstart_sup+nsup-1
2286 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2288 if (wfrag_(i).gt.0.0d0) then
2289 do j=ifrag_(1,i),ifrag_(2,i)-1
2290 do k=j+1,ifrag_(2,i)
2291 c write (iout,*) "j",j," k",k
2293 if (constr_dist.eq.1) then
2298 forcon(nhpb)=wfrag_(i)
2299 else if (constr_dist.eq.2) then
2300 if (ddjk.le.dist_cut) then
2305 forcon(nhpb)=wfrag_(i)
2312 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2315 if (.not.out1file .or. me.eq.king)
2316 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2317 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2319 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2320 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2327 if (wpair_(i).gt.0.0d0) then
2335 do j=ifrag_(1,ii),ifrag_(2,ii)
2336 do k=ifrag_(1,jj),ifrag_(2,jj)
2340 forcon(nhpb)=wpair_(i)
2341 dhpb(nhpb)=dist(j,k)
2343 if (.not.out1file .or. me.eq.king)
2344 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2345 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2347 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2348 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2355 if (constr_dist.eq.11) then
2356 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2357 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2358 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2360 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2361 & ibecarb(i),forcon(nhpb+1)
2363 if (forcon(nhpb+1).gt.0.0d0) then
2365 if (ibecarb(i).gt.0) then
2366 ihpb(i)=ihpb(i)+nres
2367 jhpb(i)=jhpb(i)+nres
2369 if (dhpb(nhpb).eq.0.0d0)
2370 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2372 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2373 C if (forcon(nhpb+1).gt.0.0d0) then
2375 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2377 if (.not.out1file .or. me.eq.king)
2378 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2379 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2381 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2382 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2389 c-------------------------------------------------------------------------------
2391 subroutine flush(iu)
2396 subroutine flush(iu)
2401 c------------------------------------------------------------------------------
2402 subroutine copy_to_tmp(source)
2403 include "DIMENSIONS"
2404 include "COMMON.IOUNITS"
2405 character*(*) source
2406 character* 256 tmpfile
2410 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2411 inquire(file=tmpfile,exist=ex)
2413 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2414 & " to temporary directory..."
2415 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2416 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2420 c------------------------------------------------------------------------------
2421 subroutine move_from_tmp(source)
2422 include "DIMENSIONS"
2423 include "COMMON.IOUNITS"
2424 character*(*) source
2427 write (*,*) "Moving ",source(:ilen(source)),
2428 & " from temporary directory to working directory"
2429 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2430 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2433 c------------------------------------------------------------------------------
2434 subroutine random_init(seed)
2436 C Initialize random number generator
2438 implicit real*8 (a-h,o-z)
2439 include 'DIMENSIONS'
2442 logical OKRandom, prng_restart
2444 integer iseed_array(4)
2446 include 'COMMON.IOUNITS'
2447 include 'COMMON.TIME1'
2448 include 'COMMON.THREAD'
2449 include 'COMMON.SBRIDGE'
2450 include 'COMMON.CONTROL'
2451 include 'COMMON.MCM'
2452 include 'COMMON.MAP'
2453 include 'COMMON.HEADER'
2454 include 'COMMON.CSA'
2455 include 'COMMON.CHAIN'
2456 include 'COMMON.MUCA'
2458 include 'COMMON.FFIELD'
2459 include 'COMMON.SETUP'
2460 iseed=-dint(dabs(seed))
2461 if (iseed.eq.0) then
2462 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2463 & 'Random seed undefined. The program will stop.'
2464 write (*,'(/80(1h*)/20x,a/80(1h*))')
2465 & 'Random seed undefined. The program will stop.'
2467 call mpi_finalize(mpi_comm_world,ierr)
2469 stop 'Bad random seed.'
2472 if (fg_rank.eq.0) then
2476 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2477 OKRandom = prng_restart(me,iseed)
2480 tmp=65536.0d0**(4-i)
2481 iseed_array(i) = dint(seed/tmp)
2482 seed=seed-iseed_array(i)*tmp
2485 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2486 & (iseed_array(i),i=1,4)
2487 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2488 & (iseed_array(i),i=1,4)
2489 OKRandom = prng_restart(me,iseed_array)
2492 c r1 = prng_next(me)
2493 r1=ran_number(0.0D0,1.0D0)
2495 & write (iout,*) 'ran_num',r1
2496 if (r1.lt.0.0d0) OKRandom=.false.
2498 if (.not.OKRandom) then
2499 write (iout,*) 'PRNG IS NOT WORKING!!!'
2500 print *,'PRNG IS NOT WORKING!!!'
2503 call mpi_abort(mpi_comm_world,error_msg,ierr)
2506 write (iout,*) 'too many processors for parallel prng'
2507 write (*,*) 'too many processors for parallel prng'
2515 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)