2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 csa if (modecalc.eq.8) then
33 csa inquire (file="fort.40",exist=file_exist)
34 csa if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
49 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
54 c print *,"Processor",myrank," leaves READRTNS"
57 C-------------------------------------------------------------------------------
58 subroutine read_control
62 implicit real*8 (a-h,o-z)
66 logical OKRandom, prng_restart
69 include 'COMMON.IOUNITS'
70 include 'COMMON.TIME1'
71 include 'COMMON.THREAD'
72 include 'COMMON.SBRIDGE'
73 include 'COMMON.CONTROL'
76 include 'COMMON.HEADER'
77 csa include 'COMMON.CSA'
78 include 'COMMON.CHAIN'
81 include 'COMMON.FFIELD'
82 include 'COMMON.SETUP'
83 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
84 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
86 character*320 controlcard
91 read (INP,'(a)') titel
92 call card_concat(controlcard)
93 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
94 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
95 call reada(controlcard,'SEED',seed,0.0D0)
96 call random_init(seed)
97 C Set up the time limit (caution! The time must be input in minutes!)
98 read_cart=index(controlcard,'READ_CART').gt.0
99 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
100 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
101 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
102 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
103 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
104 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
105 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
106 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
107 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
108 call reada(controlcard,'DRMS',drms,0.1D0)
109 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
110 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
111 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
112 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
113 write (iout,'(a,f10.1)')'DRMS = ',drms
114 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
115 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
117 call readi(controlcard,'NZ_START',nz_start,0)
118 call readi(controlcard,'NZ_END',nz_end,0)
119 call readi(controlcard,'IZ_SC',iz_sc,0)
121 safety = 60.0d0*safety
124 call reada(controlcard,"T_BATH",t_bath,300.0d0)
125 minim=(index(controlcard,'MINIMIZE').gt.0)
126 dccart=(index(controlcard,'CART').gt.0)
127 overlapsc=(index(controlcard,'OVERLAP').gt.0)
128 overlapsc=.not.overlapsc
129 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
130 searchsc=.not.searchsc
131 sideadd=(index(controlcard,'SIDEADD').gt.0)
132 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
133 outpdb=(index(controlcard,'PDBOUT').gt.0)
134 outx=(index(controlcard,'XOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 if (index(controlcard,'CHECKGRAD').gt.0) then
155 if (index(controlcard,'CART').gt.0) then
157 elseif (index(controlcard,'CARINT').gt.0) then
162 elseif (index(controlcard,'THREAD').gt.0) then
164 call readi(controlcard,'THREAD',nthread,0)
165 if (nthread.gt.0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
169 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
170 stop 'Error termination in Read_Control.'
172 else if (index(controlcard,'MCMA').gt.0) then
174 else if (index(controlcard,'MCEE').gt.0) then
176 else if (index(controlcard,'MULTCONF').gt.0) then
178 else if (index(controlcard,'MAP').gt.0) then
180 call readi(controlcard,'MAP',nmap,0)
181 else if (index(controlcard,'CSA').gt.0) then
182 write(*,*) "CSA not supported in this version"
185 crc else if (index(controlcard,'ZSCORE').gt.0) then
187 crc ZSCORE is rm from UNRES, modecalc=9 is available
190 cfcm else if (index(controlcard,'MCMF').gt.0) then
192 else if (index(controlcard,'SOFTREG').gt.0) then
194 else if (index(controlcard,'CHECK_BOND').gt.0) then
196 else if (index(controlcard,'TEST').gt.0) then
198 else if (index(controlcard,'MD').gt.0) then
200 else if (index(controlcard,'RE ').gt.0) then
204 lmuca=index(controlcard,'MUCA').gt.0
205 call readi(controlcard,'MUCADYN',mucadyn,0)
206 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
207 if (lmuca .and. (me.eq.king .or. .not.out1file ))
209 write (iout,*) 'MUCADYN=',mucadyn
210 write (iout,*) 'MUCASMOOTH=',muca_smooth
213 iscode=index(controlcard,'ONE_LETTER')
214 indphi=index(controlcard,'PHI')
215 indback=index(controlcard,'BACK')
216 iranconf=index(controlcard,'RAND_CONF')
217 i2ndstr=index(controlcard,'USE_SEC_PRED')
218 gradout=index(controlcard,'GRADOUT').gt.0
219 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
221 if(me.eq.king.or..not.out1file)
222 & write (iout,'(2a)') diagmeth(kdiag),
223 & ' routine used to diagonalize matrices.'
226 c--------------------------------------------------------------------------
227 subroutine read_REMDpar
231 implicit real*8 (a-h,o-z)
233 include 'COMMON.IOUNITS'
234 include 'COMMON.TIME1'
237 include 'COMMON.LANGEVIN'
239 include 'COMMON.LANGEVIN.lang0'
241 include 'COMMON.INTERACT'
242 include 'COMMON.NAMES'
244 include 'COMMON.REMD'
245 include 'COMMON.CONTROL'
246 include 'COMMON.SETUP'
248 character*320 controlcard
249 character*3200 controlcard1
250 integer iremd_m_total
252 if(me.eq.king.or..not.out1file)
253 & write (iout,*) "REMD setup"
255 call card_concat(controlcard)
256 call readi(controlcard,"NREP",nrep,3)
257 call readi(controlcard,"NSTEX",nstex,1000)
258 call reada(controlcard,"RETMIN",retmin,10.0d0)
259 call reada(controlcard,"RETMAX",retmax,1000.0d0)
260 mremdsync=(index(controlcard,'SYNC').gt.0)
261 call readi(controlcard,"NSYN",i_sync_step,100)
262 restart1file=(index(controlcard,'REST1FILE').gt.0)
263 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
264 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
265 if(max_cache_traj_use.gt.max_cache_traj)
266 & max_cache_traj_use=max_cache_traj
267 if(me.eq.king.or..not.out1file) then
268 cd if (traj1file) then
269 crc caching is in testing - NTWX is not ignored
270 cd write (iout,*) "NTWX value is ignored"
271 cd write (iout,*) " trajectory is stored to one file by master"
272 cd write (iout,*) " before exchange at NSTEX intervals"
274 write (iout,*) "NREP= ",nrep
275 write (iout,*) "NSTEX= ",nstex
276 write (iout,*) "SYNC= ",mremdsync
277 write (iout,*) "NSYN= ",i_sync_step
278 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
281 t_exchange_only=(index(controlcard,'TONLY').gt.0)
282 call readi(controlcard,"HREMD",hremd,0)
283 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
284 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
286 if(usampl.and.hremd.gt.0) then
288 & "========== ERROR: USAMPL and HREMD cannot be used together"
290 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
297 if (index(controlcard,'TLIST').gt.0) then
299 call card_concat(controlcard1)
300 read(controlcard1,*) (remd_t(i),i=1,nrep)
301 if(me.eq.king.or..not.out1file)
302 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
305 if (index(controlcard,'MLIST').gt.0) then
307 call card_concat(controlcard1)
308 read(controlcard1,*) (remd_m(i),i=1,nrep)
309 if(me.eq.king.or..not.out1file) then
310 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
313 iremd_m_total=iremd_m_total+remd_m(i)
316 write (iout,*) 'Total number of replicas ',
317 & iremd_m_total*hremd
319 write (iout,*) 'Total number of replicas ',iremd_m_total
323 if(me.eq.king.or..not.out1file)
324 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
327 c--------------------------------------------------------------------------
328 subroutine read_MDpar
332 implicit real*8 (a-h,o-z)
334 include 'COMMON.IOUNITS'
335 include 'COMMON.TIME1'
338 include 'COMMON.LANGEVIN'
340 include 'COMMON.LANGEVIN.lang0'
342 include 'COMMON.INTERACT'
343 include 'COMMON.NAMES'
345 include 'COMMON.SETUP'
346 include 'COMMON.CONTROL'
347 include 'COMMON.SPLITELE'
349 character*320 controlcard
351 call card_concat(controlcard)
352 call readi(controlcard,"NSTEP",n_timestep,1000000)
353 call readi(controlcard,"NTWE",ntwe,100)
354 call readi(controlcard,"NTWX",ntwx,1000)
355 call reada(controlcard,"DT",d_time,1.0d-1)
356 call reada(controlcard,"DVMAX",dvmax,2.0d1)
357 call reada(controlcard,"DAMAX",damax,1.0d1)
358 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
359 call readi(controlcard,"LANG",lang,0)
360 RESPA = index(controlcard,"RESPA") .gt. 0
361 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
362 ntime_split0=ntime_split
363 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
364 ntime_split0=ntime_split
365 call reada(controlcard,"R_CUT",r_cut,2.0d0)
366 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
367 rest = index(controlcard,"REST").gt.0
368 tbf = index(controlcard,"TBF").gt.0
369 call readi(controlcard,"HMC",hmc,0)
370 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
371 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
372 tnh = index(controlcard,"NOSEHOOVER96").gt.0
373 if (RESPA.and.tnh)then
374 xiresp = index(controlcard,"XIRESP").gt.0
376 call reada(controlcard,"Q_NP",Q_np,0.1d0)
377 usampl = index(controlcard,"USAMPL").gt.0
378 mdpdb = index(controlcard,"MDPDB").gt.0
379 call reada(controlcard,"T_BATH",t_bath,300.0d0)
380 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
381 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
382 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
383 if (count_reset_moment.eq.0) count_reset_moment=1000000000
384 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
385 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
386 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
387 if (count_reset_vel.eq.0) count_reset_vel=1000000000
388 large = index(controlcard,"LARGE").gt.0
389 print_compon = index(controlcard,"PRINT_COMPON").gt.0
390 rattle = index(controlcard,"RATTLE").gt.0
391 preminim = index(controlcard,"PREMINIM").gt.0
393 dccart=(index(controlcard,'CART').gt.0)
396 c if performing umbrella sampling, fragments constrained are read from the fragment file
402 if(me.eq.king.or..not.out1file) then
404 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
406 write (iout,'(a)') "The units are:"
407 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
408 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
409 & " acceleration: angstrom/(48.9 fs)**2"
410 write (iout,'(a)') "energy: kcal/mol, temperature: K"
412 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
413 write (iout,'(a60,f10.5,a)')
414 & "Initial time step of numerical integration:",d_time,
416 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
418 write (iout,'(2a,i4,a)')
419 & "A-MTS algorithm used; initial time step for fast-varying",
420 & " short-range forces split into",ntime_split," steps."
421 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
422 & r_cut," lambda",rlamb
424 write (iout,'(2a,f10.5)')
425 & "Maximum acceleration threshold to reduce the time step",
426 & "/increase split number:",damax
427 write (iout,'(2a,f10.5)')
428 & "Maximum predicted energy drift to reduce the timestep",
429 & "/increase split number:",edriftmax
430 write (iout,'(a60,f10.5)')
431 & "Maximum velocity threshold to reduce velocities:",dvmax
432 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
433 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
434 if (rattle) write (iout,'(a60)')
435 & "Rattle algorithm used to constrain the virtual bonds"
436 if (preminim .or. iranconf.gt.0) then
438 & "Initial structure will be energy-minimized"
443 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
444 call reada(controlcard,"RWAT",rwat,1.4d0)
445 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
446 surfarea=index(controlcard,"SURFAREA").gt.0
447 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
448 if(me.eq.king.or..not.out1file)then
449 write (iout,'(/a,$)') "Langevin dynamics calculation"
452 & " with direct integration of Langevin equations"
453 else if (lang.eq.2) then
454 write (iout,'(a/)') " with TINKER stochasic MD integrator"
455 else if (lang.eq.3) then
456 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
457 else if (lang.eq.4) then
458 write (iout,'(a/)') " in overdamped mode"
460 write (iout,'(//a,i5)')
461 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
464 write (iout,'(a60,f10.5)') "Temperature:",t_bath
465 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
466 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
467 write (iout,'(a60,f10.5)')
468 & "Scaling factor of the friction forces:",scal_fric
469 if (surfarea) write (iout,'(2a,i10,a)')
470 & "Friction coefficients will be scaled by solvent-accessible",
471 & " surface area every",reset_fricmat," steps."
473 c Calculate friction coefficients and bounds of stochastic forces
474 eta=6*pi*cPoise*etawat
475 if(me.eq.king.or..not.out1file)
476 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
478 gamp=scal_fric*(pstok+rwat)*eta
479 stdfp=dsqrt(2*Rb*t_bath/d_time)
481 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
482 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(/2a/)')
486 & "Radii of site types and friction coefficients and std's of",
487 & " stochastic forces of fully exposed sites"
488 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
490 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
491 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
495 if(me.eq.king.or..not.out1file)then
496 write (iout,'(a)') "Berendsen bath calculation"
497 write (iout,'(a60,f10.5)') "Temperature:",t_bath
498 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
500 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
501 & count_reset_moment," steps"
503 & write (iout,'(a,i10,a)')
504 & "Velocities will be reset at random every",count_reset_vel,
507 else if (tnp .or. tnp1 .or. tnh) then
508 if (tnp .or. tnp1) then
509 write (iout,'(a)') "Nose-Poincare bath calculation"
510 if (tnp) write (iout,'(a)')
511 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
512 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
514 write (iout,'(a)') "Nose-Hoover bath calculation"
515 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
525 WDTI(i) = 1.0*d_time/nresn
532 write (iout,'(a)') "NVT-XI-RESPA algorithm"
534 write (iout,'(a)') "NVT-XO-RESPA algorithm"
537 WDTIi(i) = 1.0*d_time/nresn/ntime_split
545 write (iout,'(a60,f10.5)') "Temperature:",t_bath
546 write (iout,'(a60,f10.5)') "Q =",Q_np
548 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
549 & count_reset_moment," steps"
551 & write (iout,'(a,i10,a)')
552 & "Velocities will be reset at random every",count_reset_vel,
555 else if (hmc.gt.0) then
556 write (iout,'(a)') "Hybrid Monte Carlo calculation"
557 write (iout,'(a60,f10.5)') "Temperature:",t_bath
558 write (iout,'(a60,i10)')
559 & "Number of MD steps between Metropolis tests:",hmc
562 if(me.eq.king.or..not.out1file)
563 & write (iout,'(a31)') "Microcanonical mode calculation"
565 if(me.eq.king.or..not.out1file)then
566 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
568 write(iout,*) "MD running with constraints."
569 write(iout,*) "Equilibration time ", eq_time, " mtus."
570 write(iout,*) "Constraining ", nfrag," fragments."
571 write(iout,*) "Length of each fragment, weight and q0:"
573 write (iout,*) "Set of restraints #",iset
575 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
576 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
578 write(iout,*) "constraints between ", npair, "fragments."
579 write(iout,*) "constraint pairs, weights and q0:"
581 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
582 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
584 write(iout,*) "angle constraints within ", nfrag_back,
585 & "backbone fragments."
586 write(iout,*) "fragment, weights:"
588 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
589 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
590 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
593 iset=mod(kolor,nset)+1
596 if(me.eq.king.or..not.out1file)
597 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
600 c------------------------------------------------------------------------------
603 C Read molecular data.
605 implicit real*8 (a-h,o-z)
611 include 'COMMON.IOUNITS'
614 include 'COMMON.INTERACT'
615 include 'COMMON.LOCAL'
616 include 'COMMON.NAMES'
617 include 'COMMON.CHAIN'
618 include 'COMMON.FFIELD'
619 include 'COMMON.SBRIDGE'
620 include 'COMMON.HEADER'
621 include 'COMMON.CONTROL'
622 include 'COMMON.DBASE'
623 include 'COMMON.THREAD'
624 include 'COMMON.CONTACTS'
625 include 'COMMON.TORCNSTR'
626 include 'COMMON.TIME1'
627 include 'COMMON.BOUNDS'
629 include 'COMMON.REMD'
630 include 'COMMON.SETUP'
631 character*4 sequence(maxres)
633 double precision x(maxvar)
634 character*256 pdbfile
635 character*320 weightcard
636 character*80 weightcard_t,ucase
637 dimension itype_pdb(maxres)
638 common /pizda/ itype_pdb
639 logical seq_comp,fail
640 double precision energia(0:n_ene)
647 C Read weights of the subsequent energy terms.
660 if(me.eq.king.or..not.out1file) then
661 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
662 write (iout,*) 'Current weights for processor ',
663 & me,' set ',i2set(me)
667 call card_concat(weightcard)
668 call reada(weightcard,'WLONG',wlong,1.0D0)
669 call reada(weightcard,'WSC',wsc,wlong)
670 call reada(weightcard,'WSCP',wscp,wlong)
671 call reada(weightcard,'WELEC',welec,1.0D0)
672 call reada(weightcard,'WVDWPP',wvdwpp,welec)
673 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
674 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
675 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
676 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
677 call reada(weightcard,'WTURN3',wturn3,1.0D0)
678 call reada(weightcard,'WTURN4',wturn4,1.0D0)
679 call reada(weightcard,'WTURN6',wturn6,1.0D0)
680 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
681 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
682 call reada(weightcard,'WBOND',wbond,1.0D0)
683 call reada(weightcard,'WTOR',wtor,1.0D0)
684 call reada(weightcard,'WTORD',wtor_d,1.0D0)
685 call reada(weightcard,'WANG',wang,1.0D0)
686 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
687 call reada(weightcard,'SCAL14',scal14,0.4D0)
688 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
689 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
690 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
691 call reada(weightcard,'TEMP0',temp0,300.0d0)
692 if (index(weightcard,'SOFT').gt.0) ipot=6
693 C 12/1/95 Added weight for the multi-body term WCORR
694 call reada(weightcard,'WCORRH',wcorr,1.0D0)
695 if (wcorr4.gt.0.0d0) wcorr=wcorr4
703 hweights(i,7)=wel_loc
706 hweights(i,10)=wturn6
708 hweights(i,12)=wscloc
710 hweights(i,14)=wtor_d
711 hweights(i,15)=wstrain
712 hweights(i,16)=wvdwpp
714 hweights(i,18)=scal14
715 hweights(i,21)=wsccor
720 weights(i)=hweights(i2set(me),i)
744 call card_concat(weightcard)
745 call reada(weightcard,'WLONG',wlong,1.0D0)
746 call reada(weightcard,'WSC',wsc,wlong)
747 call reada(weightcard,'WSCP',wscp,wlong)
748 call reada(weightcard,'WELEC',welec,1.0D0)
749 call reada(weightcard,'WVDWPP',wvdwpp,welec)
750 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
751 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
752 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
753 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
754 call reada(weightcard,'WTURN3',wturn3,1.0D0)
755 call reada(weightcard,'WTURN4',wturn4,1.0D0)
756 call reada(weightcard,'WTURN6',wturn6,1.0D0)
757 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
758 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
759 call reada(weightcard,'WBOND',wbond,1.0D0)
760 call reada(weightcard,'WTOR',wtor,1.0D0)
761 call reada(weightcard,'WTORD',wtor_d,1.0D0)
762 call reada(weightcard,'WANG',wang,1.0D0)
763 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
764 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
765 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
766 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
767 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
768 call reada(weightcard,'SCAL14',scal14,0.4D0)
769 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
770 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
771 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
772 call reada(weightcard,'TEMP0',temp0,300.0d0)
773 if (index(weightcard,'SOFT').gt.0) ipot=6
774 C 12/1/95 Added weight for the multi-body term WCORR
775 call reada(weightcard,'WCORRH',wcorr,1.0D0)
776 if (wcorr4.gt.0.0d0) wcorr=wcorr4
797 weights(25)=wdfa_dist
800 weights(28)=wdfa_beta
802 if(me.eq.king.or..not.out1file)
803 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
804 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
806 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
807 10 format (/'Energy-term weights (unscaled):'//
808 & 'WSCC= ',f10.6,' (SC-SC)'/
809 & 'WSCP= ',f10.6,' (SC-p)'/
810 & 'WELEC= ',f10.6,' (p-p electr)'/
811 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
812 & 'WBOND= ',f10.6,' (stretching)'/
813 & 'WANG= ',f10.6,' (bending)'/
814 & 'WSCLOC= ',f10.6,' (SC local)'/
815 & 'WTOR= ',f10.6,' (torsional)'/
816 & 'WTORD= ',f10.6,' (double torsional)'/
817 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
818 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
819 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
820 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
821 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
822 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
823 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
824 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
825 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
826 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
827 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
828 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
829 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
830 if(me.eq.king.or..not.out1file)then
831 if (wcorr4.gt.0.0d0) then
832 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
833 & 'between contact pairs of peptide groups'
834 write (iout,'(2(a,f5.3/))')
835 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
836 & 'Range of quenching the correlation terms:',2*delt_corr
837 else if (wcorr.gt.0.0d0) then
838 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
839 & 'between contact pairs of peptide groups'
841 write (iout,'(a,f8.3)')
842 & 'Scaling factor of 1,4 SC-p interactions:',scal14
843 write (iout,'(a,f8.3)')
844 & 'General scaling factor of SC-p interactions:',scalscp
846 r0_corr=cutoff_corr-delt_corr
848 aad(i,1)=scalscp*aad(i,1)
849 aad(i,2)=scalscp*aad(i,2)
850 bad(i,1)=scalscp*bad(i,1)
851 bad(i,2)=scalscp*bad(i,2)
853 call rescale_weights(t_bath)
854 if(me.eq.king.or..not.out1file)
855 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
856 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
858 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
859 22 format (/'Energy-term weights (scaled):'//
860 & 'WSCC= ',f10.6,' (SC-SC)'/
861 & 'WSCP= ',f10.6,' (SC-p)'/
862 & 'WELEC= ',f10.6,' (p-p electr)'/
863 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
864 & 'WBOND= ',f10.6,' (stretching)'/
865 & 'WANG= ',f10.6,' (bending)'/
866 & 'WSCLOC= ',f10.6,' (SC local)'/
867 & 'WTOR= ',f10.6,' (torsional)'/
868 & 'WTORD= ',f10.6,' (double torsional)'/
869 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
870 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
871 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
872 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
873 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
874 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
875 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
876 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
877 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
878 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
879 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
880 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
881 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
882 if(me.eq.king.or..not.out1file)
883 & write (iout,*) "Reference temperature for weights calculation:",
885 call reada(weightcard,"D0CM",d0cm,3.78d0)
886 call reada(weightcard,"AKCM",akcm,15.1d0)
887 call reada(weightcard,"AKTH",akth,11.0d0)
888 call reada(weightcard,"AKCT",akct,12.0d0)
889 call reada(weightcard,"V1SS",v1ss,-1.08d0)
890 call reada(weightcard,"V2SS",v2ss,7.61d0)
891 call reada(weightcard,"V3SS",v3ss,13.7d0)
892 call reada(weightcard,"EBR",ebr,-5.50D0)
893 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
895 dyn_ss_mask(i)=.false.
899 dyn_ssbond_ij(i,j)=1.0d300
902 call reada(weightcard,"HT",Ht,0.0D0)
904 ss_depth=ebr/wsc-0.25*eps(1,1)
905 Ht=Ht/wsc-0.25*eps(1,1)
906 akcm=akcm*wstrain/wsc
907 akth=akth*wstrain/wsc
908 akct=akct*wstrain/wsc
909 v1ss=v1ss*wstrain/wsc
910 v2ss=v2ss*wstrain/wsc
911 v3ss=v3ss*wstrain/wsc
913 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
916 if(me.eq.king.or..not.out1file) then
917 write (iout,*) "Parameters of the SS-bond potential:"
918 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
920 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
921 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
922 write (iout,*)" HT",Ht
923 print *,'indpdb=',indpdb,' pdbref=',pdbref
925 if (indpdb.gt.0 .or. pdbref) then
926 read(inp,'(a)') pdbfile
927 if(me.eq.king.or..not.out1file)
928 & write (iout,'(2a)') 'PDB data will be read from file ',
929 & pdbfile(:ilen(pdbfile))
930 open(ipdbin,file=pdbfile,status='old',err=33)
932 33 write (iout,'(a)') 'Error opening PDB file.'
935 c print *,'Begin reading pdb data'
944 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
945 & (crefjlee(j,i+nres),j=1,3)
948 c print *,'Finished reading pdb data'
949 if(me.eq.king.or..not.out1file)
950 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
951 & ' nstart_sup=',nstart_sup
953 itype_pdb(i)=itype(i)
957 nct=nstart_sup+nsup-1
958 call contact(.false.,ncont_ref,icont_ref,co)
961 C Following 2 lines for diagnostics; comment out if not needed
962 write (iout,*) "Before sideadd"
964 if(me.eq.king.or..not.out1file)
965 & write(iout,*)'Adding sidechains'
972 do while (fail.and.nsi.le.maxsi)
973 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
976 if(fail) write(iout,*)'Adding sidechain failed for res ',
977 & i,' after ',nsi,' trials'
980 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
983 C Following 2 lines for diagnostics; comment out if not needed
984 c write (iout,*) "After sideadd"
987 if (indpdb.eq.0) then
988 C Read sequence if not taken from the pdb file.
990 c print *,'nres=',nres
991 if (iscode.gt.0) then
992 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
994 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
996 C Convert sequence to numeric code
998 itype(i)=rescode(i,sequence(i),iscode)
1000 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
1002 & "Glycine is the first full residue, initial dummy deleted"
1008 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
1010 & "Glycine is the last full residue, terminal dummy deleted"
1014 C Assign initial virtual bond lengths
1020 vbld(i+nres)=dsc(itype(i))
1021 vbld_inv(i+nres)=dsc_inv(itype(i))
1022 c write (iout,*) "i",i," itype",itype(i),
1023 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1027 c print '(20i4)',(itype(i),i=1,nres)
1030 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1032 if (itype(i).eq.21) then
1036 else if (itype(i+1).ne.20) then
1038 else if (itype(i).ne.20) then
1045 if(me.eq.king.or..not.out1file)then
1046 write (iout,*) "ITEL"
1048 write (iout,*) i,itype(i),itel(i)
1050 print *,'Call Read_Bridge.'
1053 C 8/13/98 Set limits to generating the dihedral angles
1058 read (inp,*) ndih_constr
1059 if (ndih_constr.gt.0) then
1061 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1062 if(me.eq.king.or..not.out1file)then
1064 & 'There are',ndih_constr,' constraints on phi angles.'
1066 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1070 phi0(i)=deg2rad*phi0(i)
1071 drange(i)=deg2rad*drange(i)
1073 if(me.eq.king.or..not.out1file)
1074 & write (iout,*) 'FTORS',ftors
1077 phibound(1,ii) = phi0(i)-drange(i)
1078 phibound(2,ii) = phi0(i)+drange(i)
1083 if (me.eq.king) then
1085 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1087 write (iout,'(a3,i5,2f10.1)')
1088 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1094 cd print *,'NNT=',NNT,' NCT=',NCT
1095 if (itype(1).eq.21) nnt=2
1096 if (itype(nres).eq.21) nct=nct-1
1098 C Bartek:READ init_vars
1099 C Initialize variables!
1100 C Juyong:READ read_info
1101 C READ fragment information!!
1102 C both routines should be in dfa.F file!!
1105 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1106 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1108 print*, 'init_dfa_vars finished!'
1110 print*, 'read_dfa_info finished!'
1115 if(me.eq.king.or..not.out1file)
1116 & write (iout,'(a,i3)') 'nsup=',nsup
1118 if (nsup.le.(nct-nnt+1)) then
1119 do i=0,nct-nnt+1-nsup
1120 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1126 & 'Error - sequences to be superposed do not match.'
1129 do i=0,nsup-(nct-nnt+1)
1130 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1132 nstart_sup=nstart_sup+i
1138 & 'Error - sequences to be superposed do not match.'
1141 if (nsup.eq.0) nsup=nct-nnt
1142 if (nstart_sup.eq.0) nstart_sup=nnt
1143 if (nstart_seq.eq.0) nstart_seq=nnt
1144 if(me.eq.king.or..not.out1file)
1145 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1146 & ' nstart_seq=',nstart_seq
1148 c--- Zscore rms -------
1149 if (nz_start.eq.0) nz_start=nnt
1150 if (nz_end.eq.0 .and. nsup.gt.0) then
1152 else if (nz_end.eq.0) then
1155 if(me.eq.king.or..not.out1file)then
1156 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1157 write (iout,*) 'IZ_SC=',iz_sc
1159 c----------------------
1162 if (.not.pdbref) then
1163 call read_angles(inp,*38)
1165 38 write (iout,'(a)') 'Error reading reference structure.'
1167 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1168 stop 'Error reading reference structure'
1172 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1181 call contact(.true.,ncont_ref,icont_ref,co)
1183 if(me.eq.king.or..not.out1file)
1184 & write (iout,*) 'Contact order:',co
1186 if(me.eq.king.or..not.out1file)
1187 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1190 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1192 if(me.eq.king.or..not.out1file)
1193 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1194 & icont_ref(1,i),' ',
1195 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1199 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1200 if (constr_dist.gt.0) then
1201 call read_dist_constr
1205 if (constr_homology.gt.0) then
1206 call read_constr_homology
1207 if (indpdb.gt.0 .or. pdbref) then
1210 c(j,i)=crefjlee(j,i)
1211 cref(j,i)=crefjlee(j,i)
1216 write (iout,*) "Array C"
1218 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1219 & (c(j,i+nres),j=1,3)
1221 write (iout,*) "Array Cref"
1223 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1224 & (cref(j,i+nres),j=1,3)
1227 call int_from_cart1(.false.)
1228 call sc_loc_geom(.false.)
1230 thetaref(i)=theta(i)
1235 dc(j,i)=c(j,i+1)-c(j,i)
1236 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1241 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1242 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1250 if (nhpb.gt.0) call hpb_partition
1251 c write (iout,*) "After read_dist_constr nhpb",nhpb
1253 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1254 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1255 & modecalc.ne.10) then
1256 C If input structure hasn't been supplied from the PDB file read or generate
1258 if (iranconf.eq.0 .and. .not. extconf) then
1259 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1260 & write (iout,'(a)') 'Initial geometry will be read in.'
1262 read (inp,*) time,potE,uconst,t_bath,
1263 & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1264 read(inp,'(8f10.5)',end=36,err=36)
1265 & ((c(l,k),l=1,3),k=1,nres),
1266 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1267 call int_from_cart1(.false.)
1270 dc(j,i)=c(j,i+1)-c(j,i)
1271 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1275 if (itype(i).ne.10) then
1277 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1278 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1284 call read_angles(inp,*36)
1287 36 write (iout,'(a)') 'Error reading angle file.'
1289 call mpi_finalize( MPI_COMM_WORLD,IERR )
1291 stop 'Error reading angle file.'
1293 else if (extconf) then
1294 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1295 & write (iout,'(a)') 'Extended chain initial geometry.'
1297 theta(i)=90d0*deg2rad
1300 phi(i)=180d0*deg2rad
1303 alph(i)=110d0*deg2rad
1306 omeg(i)=-120d0*deg2rad
1309 if(me.eq.king.or..not.out1file)
1310 & write (iout,'(a)') 'Random-generated initial geometry.'
1314 if (me.eq.king .or. fg_rank.eq.0 .and. (
1315 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1319 call gen_rand_conf(itmp,*30)
1321 30 write (iout,*) 'Failed to generate random conformation',
1322 & ', itrial=',itrial
1323 write (*,*) 'Processor:',me,
1324 & ' Failed to generate random conformation',
1334 write (iout,'(a,i3,a)') 'Processor:',me,
1335 & ' error in generating random conformation.'
1336 write (*,'(a,i3,a)') 'Processor:',me,
1337 & ' error in generating random conformation.'
1340 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1347 elseif (modecalc.eq.4) then
1348 read (inp,'(a)') intinname
1349 open (intin,file=intinname,status='old',err=333)
1350 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1351 & write (iout,'(a)') 'intinname',intinname
1352 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1354 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1356 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1358 stop 'Error opening angle file.'
1362 C Generate distance constraints, if the PDB structure is to be regularized.
1363 if (nthread.gt.0) then
1364 call read_threadbase
1366 write (iout,*) "READRTNS: Calling setup_var"
1369 if (me.eq.king .or. .not. out1file)
1371 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1372 write (iout,'(/a,i3,a)')
1373 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1374 write (iout,'(20i4)') (iss(i),i=1,ns)
1376 write(iout,*)"Running with dynamic disulfide-bond formation"
1378 write (iout,'(/a/)') 'Pre-formed links are:'
1384 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1385 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1391 if (ns.gt.0.and.dyn_ss) then
1395 forcon(i-nss)=forcon(i)
1402 dyn_ss_mask(iss(i))=.true.
1405 if (i2ndstr.gt.0) call secstrp2dihc
1406 c call geom_to_var(nvar,x)
1407 c call etotal(energia(0))
1408 c call enerprint(energia(0))
1409 c call briefout(0,etot)
1411 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1412 cd write (iout,'(a)') 'Variable list:'
1413 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1415 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1416 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1417 & 'Processor',myrank,': end reading molecular data.'
1421 c--------------------------------------------------------------------------
1422 logical function seq_comp(itypea,itypeb,length)
1424 integer length,itypea(length),itypeb(length)
1427 if (itypea(i).ne.itypeb(i)) then
1435 c-----------------------------------------------------------------------------
1436 subroutine read_bridge
1437 C Read information about disulfide bridges.
1438 implicit real*8 (a-h,o-z)
1439 include 'DIMENSIONS'
1443 include 'COMMON.IOUNITS'
1444 include 'COMMON.GEO'
1445 include 'COMMON.VAR'
1446 include 'COMMON.INTERACT'
1447 include 'COMMON.LOCAL'
1448 include 'COMMON.NAMES'
1449 include 'COMMON.CHAIN'
1450 include 'COMMON.FFIELD'
1451 include 'COMMON.SBRIDGE'
1452 include 'COMMON.HEADER'
1453 include 'COMMON.CONTROL'
1454 include 'COMMON.DBASE'
1455 include 'COMMON.THREAD'
1456 include 'COMMON.TIME1'
1457 include 'COMMON.SETUP'
1458 C Read bridging residues.
1459 read (inp,*) ns,(iss(i),i=1,ns)
1461 if(me.eq.king.or..not.out1file)
1462 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1463 C Check whether the specified bridging residues are cystines.
1465 if (itype(iss(i)).ne.1) then
1466 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1467 & 'Do you REALLY think that the residue ',
1468 & restyp(itype(iss(i))),i,
1469 & ' can form a disulfide bridge?!!!'
1470 write (*,'(2a,i3,a)')
1471 & 'Do you REALLY think that the residue ',
1472 & restyp(itype(iss(i))),i,
1473 & ' can form a disulfide bridge?!!!'
1475 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1480 C Read preformed bridges.
1482 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1484 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1487 C Check if the residues involved in bridges are in the specified list of
1488 C bridging residues.
1491 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1492 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1493 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1494 & ' contains residues present in other pairs.'
1495 write (*,'(a,i3,a)') 'Disulfide pair',i,
1496 & ' contains residues present in other pairs.'
1498 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1504 if (ihpb(i).eq.iss(j)) goto 10
1506 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1509 if (jhpb(i).eq.iss(j)) goto 20
1511 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1517 ihpb(i)=ihpb(i)+nres
1518 jhpb(i)=jhpb(i)+nres
1524 c----------------------------------------------------------------------------
1525 subroutine read_x(kanal,*)
1526 implicit real*8 (a-h,o-z)
1527 include 'DIMENSIONS'
1528 include 'COMMON.GEO'
1529 include 'COMMON.VAR'
1530 include 'COMMON.CHAIN'
1531 include 'COMMON.IOUNITS'
1532 include 'COMMON.CONTROL'
1533 include 'COMMON.LOCAL'
1534 include 'COMMON.INTERACT'
1535 c Read coordinates from input
1537 read(kanal,'(8f10.5)',end=10,err=10)
1538 & ((c(l,k),l=1,3),k=1,nres),
1539 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1542 c(j,2*nres)=c(j,nres)
1544 call int_from_cart1(.false.)
1547 dc(j,i)=c(j,i+1)-c(j,i)
1548 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1552 if (itype(i).ne.10) then
1554 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1555 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1563 c----------------------------------------------------------------------------
1564 subroutine read_threadbase
1565 implicit real*8 (a-h,o-z)
1566 include 'DIMENSIONS'
1567 include 'COMMON.IOUNITS'
1568 include 'COMMON.GEO'
1569 include 'COMMON.VAR'
1570 include 'COMMON.INTERACT'
1571 include 'COMMON.LOCAL'
1572 include 'COMMON.NAMES'
1573 include 'COMMON.CHAIN'
1574 include 'COMMON.FFIELD'
1575 include 'COMMON.SBRIDGE'
1576 include 'COMMON.HEADER'
1577 include 'COMMON.CONTROL'
1578 include 'COMMON.DBASE'
1579 include 'COMMON.THREAD'
1580 include 'COMMON.TIME1'
1581 C Read pattern database for threading.
1582 read (icbase,*) nseq
1584 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1585 & nres_base(2,i),nres_base(3,i)
1586 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1588 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1589 c & nres_base(2,i),nres_base(3,i)
1590 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1594 if (weidis.eq.0.0D0) weidis=0.1D0
1603 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1604 write (iout,'(a,i5)') 'nexcl: ',nexcl
1605 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1608 c------------------------------------------------------------------------------
1609 subroutine setup_var
1610 implicit real*8 (a-h,o-z)
1611 include 'DIMENSIONS'
1612 include 'COMMON.IOUNITS'
1613 include 'COMMON.GEO'
1614 include 'COMMON.VAR'
1615 include 'COMMON.INTERACT'
1616 include 'COMMON.LOCAL'
1617 include 'COMMON.NAMES'
1618 include 'COMMON.CHAIN'
1619 include 'COMMON.FFIELD'
1620 include 'COMMON.SBRIDGE'
1621 include 'COMMON.HEADER'
1622 include 'COMMON.CONTROL'
1623 include 'COMMON.DBASE'
1624 include 'COMMON.THREAD'
1625 include 'COMMON.TIME1'
1626 C Set up variable list.
1632 if (itype(i).ne.10) then
1634 ialph(i,1)=nvar+nside
1638 if (indphi.gt.0) then
1640 else if (indback.gt.0) then
1645 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1648 c----------------------------------------------------------------------------
1649 subroutine gen_dist_constr
1650 C Generate CA distance constraints.
1651 implicit real*8 (a-h,o-z)
1652 include 'DIMENSIONS'
1653 include 'COMMON.IOUNITS'
1654 include 'COMMON.GEO'
1655 include 'COMMON.VAR'
1656 include 'COMMON.INTERACT'
1657 include 'COMMON.LOCAL'
1658 include 'COMMON.NAMES'
1659 include 'COMMON.CHAIN'
1660 include 'COMMON.FFIELD'
1661 include 'COMMON.SBRIDGE'
1662 include 'COMMON.HEADER'
1663 include 'COMMON.CONTROL'
1664 include 'COMMON.DBASE'
1665 include 'COMMON.THREAD'
1666 include 'COMMON.TIME1'
1667 dimension itype_pdb(maxres)
1668 common /pizda/ itype_pdb
1670 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1671 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1672 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1674 do i=nstart_sup,nstart_sup+nsup-1
1675 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1676 cd & ' seq_pdb', restyp(itype_pdb(i))
1677 do j=i+2,nstart_sup+nsup-1
1679 ihpb(nhpb)=i+nstart_seq-nstart_sup
1680 jhpb(nhpb)=j+nstart_seq-nstart_sup
1682 dhpb(nhpb)=dist(i,j)
1685 cd write (iout,'(a)') 'Distance constraints:'
1690 cd if (ii.gt.nres) then
1695 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1696 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1697 cd & dhpb(i),forcon(i)
1701 c----------------------------------------------------------------------------
1703 implicit real*8 (a-h,o-z)
1704 include 'DIMENSIONS'
1705 include 'COMMON.MAP'
1706 include 'COMMON.IOUNITS'
1707 character*3 angid(4) /'THE','PHI','ALP','OME'/
1708 character*80 mapcard,ucase
1710 read (inp,'(a)') mapcard
1711 mapcard=ucase(mapcard)
1712 if (index(mapcard,'PHI').gt.0) then
1714 else if (index(mapcard,'THE').gt.0) then
1716 else if (index(mapcard,'ALP').gt.0) then
1718 else if (index(mapcard,'OME').gt.0) then
1721 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1722 stop 'Error - illegal variable spec in MAP card.'
1724 call readi (mapcard,'RES1',res1(imap),0)
1725 call readi (mapcard,'RES2',res2(imap),0)
1726 if (res1(imap).eq.0) then
1727 res1(imap)=res2(imap)
1728 else if (res2(imap).eq.0) then
1729 res2(imap)=res1(imap)
1731 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1733 & 'Error - illegal definition of variable group in MAP.'
1734 stop 'Error - illegal definition of variable group in MAP.'
1736 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1737 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1738 call readi(mapcard,'NSTEP',nstep(imap),0)
1739 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1741 & 'Illegal boundary and/or step size specification in MAP.'
1742 stop 'Illegal boundary and/or step size specification in MAP.'
1747 c----------------------------------------------------------------------------
1748 csa subroutine csaread
1749 csa implicit real*8 (a-h,o-z)
1750 csa include 'DIMENSIONS'
1751 csa include 'COMMON.IOUNITS'
1752 csa include 'COMMON.GEO'
1753 csa include 'COMMON.CSA'
1754 csa include 'COMMON.BANK'
1755 csa include 'COMMON.CONTROL'
1756 csa character*80 ucase
1757 csa character*620 mcmcard
1758 csa call card_concat(mcmcard)
1760 csa call readi(mcmcard,'NCONF',nconf,50)
1761 csa call readi(mcmcard,'NADD',nadd,0)
1762 csa call readi(mcmcard,'JSTART',jstart,1)
1763 csa call readi(mcmcard,'JEND',jend,1)
1764 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1765 csa call readi(mcmcard,'N0',n0,1)
1766 csa call readi(mcmcard,'N1',n1,6)
1767 csa call readi(mcmcard,'N2',n2,4)
1768 csa call readi(mcmcard,'N3',n3,0)
1769 csa call readi(mcmcard,'N4',n4,0)
1770 csa call readi(mcmcard,'N5',n5,0)
1771 csa call readi(mcmcard,'N6',n6,10)
1772 csa call readi(mcmcard,'N7',n7,0)
1773 csa call readi(mcmcard,'N8',n8,0)
1774 csa call readi(mcmcard,'N9',n9,0)
1775 csa call readi(mcmcard,'N14',n14,0)
1776 csa call readi(mcmcard,'N15',n15,0)
1777 csa call readi(mcmcard,'N16',n16,0)
1778 csa call readi(mcmcard,'N17',n17,0)
1779 csa call readi(mcmcard,'N18',n18,0)
1781 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1783 csa call readi(mcmcard,'NDIFF',ndiff,2)
1784 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1785 csa call readi(mcmcard,'IS1',is1,1)
1786 csa call readi(mcmcard,'IS2',is2,8)
1787 csa call readi(mcmcard,'NRAN0',nran0,4)
1788 csa call readi(mcmcard,'NRAN1',nran1,2)
1789 csa call readi(mcmcard,'IRR',irr,1)
1790 csa call readi(mcmcard,'NSEED',nseed,20)
1791 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1792 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1793 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1794 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1795 csa call readi(mcmcard,'ICMAX',icmax,3)
1796 csa call readi(mcmcard,'IRESTART',irestart,0)
1797 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1800 csa call reada(mcmcard,'DELE',dele,20.0d0)
1801 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1802 csa call readi(mcmcard,'IREF',iref,0)
1803 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1804 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1805 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1806 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1807 csa write (iout,*) "NCONF_IN",nconf_in
1810 c----------------------------------------------------------------------------
1811 cfmc subroutine mcmfread
1812 cfmc implicit real*8 (a-h,o-z)
1813 cfmc include 'DIMENSIONS'
1814 cfmc include 'COMMON.MCMF'
1815 cfmc include 'COMMON.IOUNITS'
1816 cfmc include 'COMMON.GEO'
1817 cfmc character*80 ucase
1818 cfmc character*620 mcmcard
1819 cfmc call card_concat(mcmcard)
1821 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1822 cfmc write(iout,*)'MAXRANT=',maxrant
1823 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1824 cfmc write(iout,*)'MAXFAM=',maxfam
1825 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1826 cfmc write(iout,*)'NNET1=',nnet1
1827 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1828 cfmc write(iout,*)'NNET2=',nnet2
1829 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1830 cfmc write(iout,*)'NNET3=',nnet3
1831 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1832 cfmc write(iout,*)'ILASTT=',ilastt
1833 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1834 cfmc write(iout,*)'MAXSTR=',maxstr
1835 cfmc maxstr_f=maxstr/maxfam
1836 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1837 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1838 cfmc write(iout,*)'NMCMF=',nmcmf
1839 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1840 cfmc write(iout,*)'IFOCUS=',ifocus
1841 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1842 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1843 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1844 cfmc write(iout,*)'INTPRT=',intprt
1845 cfmc call readi(mcmcard,'IPRT',iprt,100)
1846 cfmc write(iout,*)'IPRT=',iprt
1847 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1848 cfmc write(iout,*)'IMAXTR=',imaxtr
1849 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1850 cfmc write(iout,*)'MAXEVEN=',maxeven
1851 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1852 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1853 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1854 cfmc write(iout,*)'INIMIN=',inimin
1855 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1856 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1857 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1858 cfmc write(iout,*)'NTHREAD=',nthread
1859 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1860 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1861 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1862 cfmc write(iout,*)'MAXPERT=',maxpert
1863 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1864 cfmc write(iout,*)'IRMSD=',irmsd
1865 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1866 cfmc write(iout,*)'DENEMIN=',denemin
1867 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1868 cfmc write(iout,*)'RCUT1S=',rcut1s
1869 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1870 cfmc write(iout,*)'RCUT1E=',rcut1e
1871 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1872 cfmc write(iout,*)'RCUT2S=',rcut2s
1873 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1874 cfmc write(iout,*)'RCUT2E=',rcut2e
1875 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1876 cfmc write(iout,*)'DPERT1=',d_pert1
1877 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1878 cfmc write(iout,*)'DPERT1A=',d_pert1a
1879 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1880 cfmc write(iout,*)'DPERT2=',d_pert2
1881 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1882 cfmc write(iout,*)'DPERT2A=',d_pert2a
1883 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1884 cfmc write(iout,*)'DPERT2B=',d_pert2b
1885 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1886 cfmc write(iout,*)'DPERT2C=',d_pert2c
1887 cfmc d_pert1=deg2rad*d_pert1
1888 cfmc d_pert1a=deg2rad*d_pert1a
1889 cfmc d_pert2=deg2rad*d_pert2
1890 cfmc d_pert2a=deg2rad*d_pert2a
1891 cfmc d_pert2b=deg2rad*d_pert2b
1892 cfmc d_pert2c=deg2rad*d_pert2c
1893 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1894 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1895 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1896 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1897 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1898 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1899 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1900 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1901 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1902 cfmc write(iout,*)'RCUTINI=',rcutini
1903 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1904 cfmc write(iout,*)'GRAT=',grat
1905 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1906 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1910 c----------------------------------------------------------------------------
1912 implicit real*8 (a-h,o-z)
1913 include 'DIMENSIONS'
1914 include 'COMMON.MCM'
1915 include 'COMMON.MCE'
1916 include 'COMMON.IOUNITS'
1918 character*320 mcmcard
1919 call card_concat(mcmcard)
1920 call readi(mcmcard,'MAXACC',maxacc,100)
1921 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1922 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1923 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1924 call readi(mcmcard,'MAXREPM',maxrepm,200)
1925 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1926 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1927 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1928 call reada(mcmcard,'E_UP',e_up,5.0D0)
1929 call reada(mcmcard,'DELTE',delte,0.1D0)
1930 call readi(mcmcard,'NSWEEP',nsweep,5)
1931 call readi(mcmcard,'NSTEPH',nsteph,0)
1932 call readi(mcmcard,'NSTEPC',nstepc,0)
1933 call reada(mcmcard,'TMIN',tmin,298.0D0)
1934 call reada(mcmcard,'TMAX',tmax,298.0D0)
1935 call readi(mcmcard,'NWINDOW',nwindow,0)
1936 call readi(mcmcard,'PRINT_MC',print_mc,0)
1937 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1938 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1939 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1940 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1941 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1942 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1943 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1944 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1945 if (nwindow.gt.0) then
1946 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1948 winlen(i)=winend(i)-winstart(i)+1
1951 if (tmax.lt.tmin) tmax=tmin
1952 if (tmax.eq.tmin) then
1956 if (nstepc.gt.0 .and. nsteph.gt.0) then
1957 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1958 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1960 C Probabilities of different move types
1961 sumpro_type(0)=0.0D0
1962 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1963 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1964 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1965 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1966 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1967 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1968 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1970 print *,'i',i,' sumprotype',sumpro_type(i)
1971 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1972 print *,'i',i,' sumprotype',sumpro_type(i)
1976 c----------------------------------------------------------------------------
1977 subroutine read_minim
1978 implicit real*8 (a-h,o-z)
1979 include 'DIMENSIONS'
1980 include 'COMMON.MINIM'
1981 include 'COMMON.IOUNITS'
1983 character*320 minimcard
1984 call card_concat(minimcard)
1985 call readi(minimcard,'MAXMIN',maxmin,2000)
1986 call readi(minimcard,'MAXFUN',maxfun,5000)
1987 call readi(minimcard,'MINMIN',minmin,maxmin)
1988 call readi(minimcard,'MINFUN',minfun,maxmin)
1989 call reada(minimcard,'TOLF',tolf,1.0D-2)
1990 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1991 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1992 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1993 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1994 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1995 & 'Options in energy minimization:'
1996 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1997 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1998 & 'MinMin:',MinMin,' MinFun:',MinFun,
1999 & ' TolF:',TolF,' RTolF:',RTolF
2002 c----------------------------------------------------------------------------
2003 subroutine read_angles(kanal,*)
2004 implicit real*8 (a-h,o-z)
2005 include 'DIMENSIONS'
2006 include 'COMMON.GEO'
2007 include 'COMMON.VAR'
2008 include 'COMMON.CHAIN'
2009 include 'COMMON.IOUNITS'
2010 include 'COMMON.CONTROL'
2011 c Read angles from input
2013 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2014 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2015 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2016 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2019 c 9/7/01 avoid 180 deg valence angle
2020 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2022 theta(i)=deg2rad*theta(i)
2023 phi(i)=deg2rad*phi(i)
2024 alph(i)=deg2rad*alph(i)
2025 omeg(i)=deg2rad*omeg(i)
2030 c----------------------------------------------------------------------------
2031 subroutine reada(rekord,lancuch,wartosc,default)
2033 character*(*) rekord,lancuch
2034 double precision wartosc,default
2037 iread=index(rekord,lancuch)
2038 if (iread.eq.0) then
2042 iread=iread+ilen(lancuch)+1
2043 read (rekord(iread:),*,err=10,end=10) wartosc
2048 c----------------------------------------------------------------------------
2049 subroutine readi(rekord,lancuch,wartosc,default)
2051 character*(*) rekord,lancuch
2052 integer wartosc,default
2055 iread=index(rekord,lancuch)
2056 if (iread.eq.0) then
2060 iread=iread+ilen(lancuch)+1
2061 read (rekord(iread:),*,err=10,end=10) wartosc
2066 c----------------------------------------------------------------------------
2067 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2070 integer tablica(dim),default
2071 character*(*) rekord,lancuch
2078 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2079 if (iread.eq.0) return
2080 iread=iread+ilen(lancuch)+1
2081 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2084 c----------------------------------------------------------------------------
2085 subroutine multreada(rekord,lancuch,tablica,dim,default)
2088 double precision tablica(dim),default
2089 character*(*) rekord,lancuch
2096 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2097 if (iread.eq.0) return
2098 iread=iread+ilen(lancuch)+1
2099 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2102 c----------------------------------------------------------------------------
2103 subroutine openunits
2104 implicit real*8 (a-h,o-z)
2105 include 'DIMENSIONS'
2108 character*16 form,nodename
2111 include 'COMMON.SETUP'
2112 include 'COMMON.IOUNITS'
2114 include 'COMMON.CONTROL'
2115 integer lenpre,lenpot,ilen,lentmp
2117 character*3 out1file_text,ucase
2120 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2121 call getenv_loc("PREFIX",prefix)
2123 call getenv_loc("POT",pot)
2124 call getenv_loc("DIRTMP",tmpdir)
2125 call getenv_loc("CURDIR",curdir)
2126 call getenv_loc("OUT1FILE",out1file_text)
2127 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2128 out1file_text=ucase(out1file_text)
2129 if (out1file_text(1:1).eq."Y") then
2132 out1file=fg_rank.gt.0
2137 if (lentmp.gt.0) then
2138 write (*,'(80(1h!))')
2139 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2140 write (*,'(80(1h!))')
2141 write (*,*)"All output files will be on node /tmp directory."
2143 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2144 if (me.eq.king) then
2145 write (*,*) "The master node is ",nodename
2146 else if (fg_rank.eq.0) then
2147 write (*,*) "I am the CG slave node ",nodename
2149 write (*,*) "I am the FG slave node ",nodename
2152 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2153 lenpre = lentmp+lenpre+1
2155 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2156 C Get the names and open the input files
2157 #if defined(WINIFL) || defined(WINPGI)
2158 open(1,file=pref_orig(:ilen(pref_orig))//
2159 & '.inp',status='old',readonly,shared)
2160 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2161 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2162 C Get parameter filenames and open the parameter files.
2163 call getenv_loc('BONDPAR',bondname)
2164 open (ibond,file=bondname,status='old',readonly,shared)
2165 call getenv_loc('THETPAR',thetname)
2166 open (ithep,file=thetname,status='old',readonly,shared)
2168 call getenv_loc('THETPARPDB',thetname_pdb)
2169 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2171 call getenv_loc('ROTPAR',rotname)
2172 open (irotam,file=rotname,status='old',readonly,shared)
2174 call getenv_loc('ROTPARPDB',rotname_pdb)
2175 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2177 call getenv_loc('TORPAR',torname)
2178 open (itorp,file=torname,status='old',readonly,shared)
2179 call getenv_loc('TORDPAR',tordname)
2180 open (itordp,file=tordname,status='old',readonly,shared)
2181 call getenv_loc('FOURIER',fouriername)
2182 open (ifourier,file=fouriername,status='old',readonly,shared)
2183 call getenv_loc('ELEPAR',elename)
2184 open (ielep,file=elename,status='old',readonly,shared)
2185 call getenv_loc('SIDEPAR',sidename)
2186 open (isidep,file=sidename,status='old',readonly,shared)
2187 #elif (defined CRAY) || (defined AIX)
2188 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2190 c print *,"Processor",myrank," opened file 1"
2191 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2192 c print *,"Processor",myrank," opened file 9"
2193 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2194 C Get parameter filenames and open the parameter files.
2195 call getenv_loc('BONDPAR',bondname)
2196 open (ibond,file=bondname,status='old',action='read')
2197 c print *,"Processor",myrank," opened file IBOND"
2198 call getenv_loc('THETPAR',thetname)
2199 open (ithep,file=thetname,status='old',action='read')
2200 c print *,"Processor",myrank," opened file ITHEP"
2202 call getenv_loc('THETPARPDB',thetname_pdb)
2203 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2205 call getenv_loc('ROTPAR',rotname)
2206 open (irotam,file=rotname,status='old',action='read')
2207 c print *,"Processor",myrank," opened file IROTAM"
2209 call getenv_loc('ROTPARPDB',rotname_pdb)
2210 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2212 call getenv_loc('TORPAR',torname)
2213 open (itorp,file=torname,status='old',action='read')
2214 c print *,"Processor",myrank," opened file ITORP"
2215 call getenv_loc('TORDPAR',tordname)
2216 open (itordp,file=tordname,status='old',action='read')
2217 c print *,"Processor",myrank," opened file ITORDP"
2218 call getenv_loc('SCCORPAR',sccorname)
2219 open (isccor,file=sccorname,status='old',action='read')
2220 c print *,"Processor",myrank," opened file ISCCOR"
2221 call getenv_loc('FOURIER',fouriername)
2222 open (ifourier,file=fouriername,status='old',action='read')
2223 c print *,"Processor",myrank," opened file IFOURIER"
2224 call getenv_loc('ELEPAR',elename)
2225 open (ielep,file=elename,status='old',action='read')
2226 c print *,"Processor",myrank," opened file IELEP"
2227 call getenv_loc('SIDEPAR',sidename)
2228 open (isidep,file=sidename,status='old',action='read')
2229 c print *,"Processor",myrank," opened file ISIDEP"
2230 c print *,"Processor",myrank," opened parameter files"
2232 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2233 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2234 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2235 C Get parameter filenames and open the parameter files.
2236 call getenv_loc('BONDPAR',bondname)
2237 open (ibond,file=bondname,status='old')
2238 call getenv_loc('THETPAR',thetname)
2239 open (ithep,file=thetname,status='old')
2241 call getenv_loc('THETPARPDB',thetname_pdb)
2242 open (ithep_pdb,file=thetname_pdb,status='old')
2244 call getenv_loc('ROTPAR',rotname)
2245 open (irotam,file=rotname,status='old')
2247 call getenv_loc('ROTPARPDB',rotname_pdb)
2248 open (irotam_pdb,file=rotname_pdb,status='old')
2250 call getenv_loc('TORPAR',torname)
2251 open (itorp,file=torname,status='old')
2252 call getenv_loc('TORDPAR',tordname)
2253 open (itordp,file=tordname,status='old')
2254 call getenv_loc('SCCORPAR',sccorname)
2255 open (isccor,file=sccorname,status='old')
2256 call getenv_loc('FOURIER',fouriername)
2257 open (ifourier,file=fouriername,status='old')
2258 call getenv_loc('ELEPAR',elename)
2259 open (ielep,file=elename,status='old')
2260 call getenv_loc('SIDEPAR',sidename)
2261 open (isidep,file=sidename,status='old')
2263 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2265 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2266 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2267 C Get parameter filenames and open the parameter files.
2268 call getenv_loc('BONDPAR',bondname)
2269 open (ibond,file=bondname,status='old',action='read')
2270 call getenv_loc('THETPAR',thetname)
2271 open (ithep,file=thetname,status='old',action='read')
2273 call getenv_loc('THETPARPDB',thetname_pdb)
2274 print *,"thetname_pdb ",thetname_pdb
2275 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2276 print *,ithep_pdb," opened"
2278 call getenv_loc('ROTPAR',rotname)
2279 open (irotam,file=rotname,status='old',action='read')
2281 call getenv_loc('ROTPARPDB',rotname_pdb)
2282 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2284 call getenv_loc('TORPAR',torname)
2285 open (itorp,file=torname,status='old',action='read')
2286 call getenv_loc('TORDPAR',tordname)
2287 open (itordp,file=tordname,status='old',action='read')
2288 call getenv_loc('SCCORPAR',sccorname)
2289 open (isccor,file=sccorname,status='old',action='read')
2290 call getenv_loc('FOURIER',fouriername)
2291 open (ifourier,file=fouriername,status='old',action='read')
2292 call getenv_loc('ELEPAR',elename)
2293 open (ielep,file=elename,status='old',action='read')
2294 call getenv_loc('SIDEPAR',sidename)
2295 open (isidep,file=sidename,status='old',action='read')
2299 C 8/9/01 In the newest version SCp interaction constants are read from a file
2300 C Use -DOLDSCP to use hard-coded constants instead.
2302 call getenv_loc('SCPPAR',scpname)
2303 #if defined(WINIFL) || defined(WINPGI)
2304 open (iscpp,file=scpname,status='old',readonly,shared)
2305 #elif (defined CRAY) || (defined AIX)
2306 open (iscpp,file=scpname,status='old',action='read')
2308 open (iscpp,file=scpname,status='old')
2310 open (iscpp,file=scpname,status='old',action='read')
2313 call getenv_loc('PATTERN',patname)
2314 #if defined(WINIFL) || defined(WINPGI)
2315 open (icbase,file=patname,status='old',readonly,shared)
2316 #elif (defined CRAY) || (defined AIX)
2317 open (icbase,file=patname,status='old',action='read')
2319 open (icbase,file=patname,status='old')
2321 open (icbase,file=patname,status='old',action='read')
2324 C Open output file only for CG processes
2325 c print *,"Processor",myrank," fg_rank",fg_rank
2326 if (fg_rank.eq.0) then
2328 if (nodes.eq.1) then
2331 npos = dlog10(dfloat(nodes-1))+1
2333 if (npos.lt.3) npos=3
2334 write (liczba,'(i1)') npos
2335 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2337 write (liczba,form) me
2338 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2339 & liczba(:ilen(liczba))
2340 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2342 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2344 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2345 & liczba(:ilen(liczba))//'.mol2'
2346 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2347 & liczba(:ilen(liczba))//'.x'
2348 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2349 & liczba(:ilen(liczba))//'.stat'
2351 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2352 & //liczba(:ilen(liczba))//'.stat')
2353 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2356 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2357 & liczba(:ilen(liczba))//'.const'
2362 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2363 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2364 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2365 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2366 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2367 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2369 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2371 rest2name=prefix(:ilen(prefix))//'.rst'
2373 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2376 #if defined(AIX) || defined(PGI)
2377 if (me.eq.king .or. .not. out1file)
2378 & open(iout,file=outname,status='unknown')
2381 if (fg_rank.gt.0) then
2382 write (liczba,'(i3.3)') myrank/nfgtasks
2383 write (ll,'(bz,i3.3)') fg_rank
2384 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2390 open(igeom,file=intname,status='unknown',position='append')
2391 open(ipdb,file=pdbname,status='unknown')
2392 open(imol2,file=mol2name,status='unknown')
2393 open(istat,file=statname,status='unknown',position='append')
2395 c1out open(iout,file=outname,status='unknown')
2398 if (me.eq.king .or. .not.out1file)
2399 & open(iout,file=outname,status='unknown')
2402 if (fg_rank.gt.0) then
2403 print "Processor",fg_rank," opening output file"
2404 write (liczba,'(i3.3)') myrank/nfgtasks
2405 write (ll,'(bz,i3.3)') fg_rank
2406 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2412 open(igeom,file=intname,status='unknown',access='append')
2413 open(ipdb,file=pdbname,status='unknown')
2414 open(imol2,file=mol2name,status='unknown')
2415 open(istat,file=statname,status='unknown',access='append')
2417 c1out open(iout,file=outname,status='unknown')
2420 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2421 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2422 csa csa_history=prefix(:lenpre)//'.CSA.history'
2423 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2424 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2425 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2426 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2427 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2428 csa csa_int=prefix(:lenpre)//'.int'
2429 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2430 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2431 csa csa_in=prefix(:lenpre)//'.CSA.in'
2432 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2435 write (iout,'(80(1h-))')
2436 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2437 write (iout,'(80(1h-))')
2438 write (iout,*) "Input file : ",
2439 & pref_orig(:ilen(pref_orig))//'.inp'
2440 write (iout,*) "Output file : ",
2441 & outname(:ilen(outname))
2443 write (iout,*) "Sidechain potential file : ",
2444 & sidename(:ilen(sidename))
2446 write (iout,*) "SCp potential file : ",
2447 & scpname(:ilen(scpname))
2449 write (iout,*) "Electrostatic potential file : ",
2450 & elename(:ilen(elename))
2451 write (iout,*) "Cumulant coefficient file : ",
2452 & fouriername(:ilen(fouriername))
2453 write (iout,*) "Torsional parameter file : ",
2454 & torname(:ilen(torname))
2455 write (iout,*) "Double torsional parameter file : ",
2456 & tordname(:ilen(tordname))
2457 write (iout,*) "SCCOR parameter file : ",
2458 & sccorname(:ilen(sccorname))
2459 write (iout,*) "Bond & inertia constant file : ",
2460 & bondname(:ilen(bondname))
2461 write (iout,*) "Bending parameter file : ",
2462 & thetname(:ilen(thetname))
2463 write (iout,*) "Rotamer parameter file : ",
2464 & rotname(:ilen(rotname))
2465 write (iout,*) "Threading database : ",
2466 & patname(:ilen(patname))
2468 &write (iout,*)" DIRTMP : ",
2470 write (iout,'(80(1h-))')
2474 c----------------------------------------------------------------------------
2475 subroutine card_concat(card)
2476 implicit real*8 (a-h,o-z)
2477 include 'DIMENSIONS'
2478 include 'COMMON.IOUNITS'
2480 character*80 karta,ucase
2482 read (inp,'(a)') karta
2485 do while (karta(80:80).eq.'&')
2486 card=card(:ilen(card)+1)//karta(:79)
2487 read (inp,'(a)') karta
2490 card=card(:ilen(card)+1)//karta
2493 c----------------------------------------------------------------------------------
2495 implicit real*8 (a-h,o-z)
2496 include 'DIMENSIONS'
2497 include 'COMMON.CHAIN'
2498 include 'COMMON.IOUNITS'
2500 include 'COMMON.CONTROL'
2501 open(irest2,file=rest2name,status='unknown')
2502 read(irest2,*) totT,EK,potE,totE,t_bath
2504 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2507 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2509 if(usampl.or.homol_nset.gt.1) then
2510 read (irest2,*) iset
2515 c---------------------------------------------------------------------------------
2516 subroutine read_fragments
2517 implicit real*8 (a-h,o-z)
2518 include 'DIMENSIONS'
2522 include 'COMMON.SETUP'
2523 include 'COMMON.CHAIN'
2524 include 'COMMON.IOUNITS'
2526 include 'COMMON.CONTROL'
2528 read(inp,*) nset,nfrag,npair,nfrag_back
2529 if(me.eq.king.or..not.out1file)
2530 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2531 & " nfrag_back",nfrag_back
2533 read(inp,*) mset(iset1)
2535 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2537 if(me.eq.king.or..not.out1file)
2538 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2539 & ifrag(2,i,iset1), qinfrag(i,iset1)
2542 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2544 if(me.eq.king.or..not.out1file)
2545 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2546 & ipair(2,i,iset1), qinpair(i,iset1)
2549 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2550 & wfrag_back(3,i,iset1),
2551 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2552 if(me.eq.king.or..not.out1file)
2553 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2554 & wfrag_back(2,i,iset1),
2555 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2556 & ifrag_back(2,i,iset1)
2561 c-------------------------------------------------------------------------------
2562 subroutine read_dist_constr
2563 implicit real*8 (a-h,o-z)
2564 include 'DIMENSIONS'
2568 include 'COMMON.SETUP'
2569 include 'COMMON.CONTROL'
2570 include 'COMMON.CHAIN'
2571 include 'COMMON.IOUNITS'
2572 include 'COMMON.SBRIDGE'
2573 integer ifrag_(2,100),ipair_(2,100)
2574 double precision wfrag_(100),wpair_(100)
2575 character*500 controlcard
2576 c write (iout,*) "Calling read_dist_constr"
2577 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2579 call card_concat(controlcard)
2580 call readi(controlcard,"NFRAG",nfrag_,0)
2581 call readi(controlcard,"NPAIR",npair_,0)
2582 call readi(controlcard,"NDIST",ndist_,0)
2583 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2585 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2586 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2587 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2588 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2589 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2590 c write (iout,*) "IFRAG"
2592 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2594 c write (iout,*) "IPAIR"
2596 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2598 if (.not.refstr .and. nfrag.gt.0) then
2600 & "ERROR: no reference structure to compute distance restraints"
2602 & "Restraints must be specified explicitly (NDIST=number)"
2605 if (nfrag.lt.2 .and. npair.gt.0) then
2606 write (iout,*) "ERROR: Less than 2 fragments specified",
2607 & " but distance restraints between pairs requested"
2612 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2613 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2614 & ifrag_(2,i)=nstart_sup+nsup-1
2615 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2617 if (wfrag_(i).gt.0.0d0) then
2618 do j=ifrag_(1,i),ifrag_(2,i)-1
2619 do k=j+1,ifrag_(2,i)
2620 c write (iout,*) "j",j," k",k
2622 if (constr_dist.eq.1) then
2627 forcon(nhpb)=wfrag_(i)
2628 else if (constr_dist.eq.2) then
2629 if (ddjk.le.dist_cut) then
2634 forcon(nhpb)=wfrag_(i)
2641 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2644 if (.not.out1file .or. me.eq.king)
2645 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2646 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2648 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2649 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2656 if (wpair_(i).gt.0.0d0) then
2664 do j=ifrag_(1,ii),ifrag_(2,ii)
2665 do k=ifrag_(1,jj),ifrag_(2,jj)
2669 forcon(nhpb)=wpair_(i)
2670 dhpb(nhpb)=dist(j,k)
2672 if (.not.out1file .or. me.eq.king)
2673 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2674 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2676 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2677 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2684 if (constr_dist.eq.11) then
2685 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2686 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2687 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2689 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2690 & ibecarb(i),forcon(nhpb+1)
2692 if (forcon(nhpb+1).gt.0.0d0) then
2694 if (ibecarb(i).gt.0) then
2695 ihpb(i)=ihpb(i)+nres
2696 jhpb(i)=jhpb(i)+nres
2698 if (dhpb(nhpb).eq.0.0d0)
2699 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2703 if (.not.out1file .or. me.eq.king) then
2706 if (constr_dist.eq.11) then
2707 write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
2708 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
2711 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2712 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2721 c-------------------------------------------------------------------------------
2723 subroutine read_constr_homology
2725 include 'DIMENSIONS'
2729 include 'COMMON.SETUP'
2730 include 'COMMON.CONTROL'
2731 include 'COMMON.CHAIN'
2732 include 'COMMON.IOUNITS'
2734 include 'COMMON.GEO'
2735 include 'COMMON.INTERACT'
2736 include 'COMMON.NAMES'
2738 c For new homol impl
2740 include 'COMMON.VAR'
2743 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2745 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2746 c & sigma_odl_temp(maxres,maxres,max_template)
2748 character*24 model_ki_dist, model_ki_angle
2749 character*500 controlcard
2750 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2751 logical lprn /.true./
2756 c FP - Nov. 2014 Temporary specifications for new vars
2758 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2759 double precision, dimension (max_template,maxres) :: rescore
2760 double precision, dimension (max_template,maxres) :: rescore2
2761 character*24 pdbfile,tpl_k_rescore
2762 c -----------------------------------------------------------------
2763 c Reading multiple PDB ref structures and calculation of retraints
2764 c not using pre-computed ones stored in files model_ki_{dist,angle}
2766 c -----------------------------------------------------------------
2769 c Alternative: reading from input
2770 call card_concat(controlcard)
2771 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2772 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2773 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2774 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2775 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2776 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2777 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2778 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2779 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2780 if(.not.read2sigma.and.start_from_model) then
2781 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2782 start_from_model=.false.
2784 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2785 if(start_from_model .and. rest) then
2786 write(iout,*) 'START_FROM_MODELS is OFF'
2787 write(iout,*) 'remove restart keyword from input'
2789 if (homol_nset.gt.1)then
2790 call card_concat(controlcard)
2791 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2792 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2793 write(iout,*) "iset homology_weight "
2795 write(iout,*) i,waga_homology(i)
2798 iset=mod(kolor,homol_nset)+1
2801 waga_homology(1)=1.0
2804 cd write (iout,*) "nnt",nnt," nct",nct
2811 write(iout,*) 'nnt=',nnt,'nct=',nct
2814 do k=1,constr_homology
2827 do k=1,constr_homology
2829 read(inp,'(a)') pdbfile
2830 c Next stament causes error upon compilation (?)
2831 c if(me.eq.king.or. .not. out1file)
2832 c write (iout,'(2a)') 'PDB data will be read from file ',
2833 c & pdbfile(:ilen(pdbfile))
2834 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2835 & pdbfile(:ilen(pdbfile))
2836 open(ipdbin,file=pdbfile,status='old',err=33)
2838 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2839 & pdbfile(:ilen(pdbfile))
2842 c print *,'Begin reading pdb data'
2844 c Files containing res sim or local scores (former containing sigmas)
2847 write(kic2,'(bz,i2.2)') k
2849 tpl_k_rescore="template"//kic2//".sco"
2852 if (read2sigma) then
2853 call readpdb_template(k)
2858 c Distance restraints
2861 C Copy the coordinates from reference coordinates (?)
2865 c write (iout,*) "c(",j,i,") =",c(j,i)
2869 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2871 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2872 open (ientin,file=tpl_k_rescore,status='old')
2873 if (nnt.gt.1) rescore(k,1)=0.0d0
2874 do irec=nnt,nct ! loop for reading res sim
2875 if (read2sigma) then
2876 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2879 idomain(k,i_tmp)=idomain_tmp
2880 rescore(k,i_tmp)=rescore_tmp
2881 rescore2(k,i_tmp)=rescore2_tmp
2882 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
2883 & i_tmp,rescore2_tmp,rescore_tmp,
2887 read (ientin,*,end=1401) rescore_tmp
2889 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2890 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2891 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2896 if (waga_dist.ne.0.0d0) then
2904 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2905 c write (iout,*) k,i,j,distal,dist2_cut
2907 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2908 & .and. distal.le.dist2_cut ) then
2914 c write (iout,*) "k",k
2915 c write (iout,*) "i",i," j",j," constr_homology",
2920 if (read2sigma) then
2923 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2925 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2926 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2927 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2929 if (odl(k,ii).le.dist_cut) then
2930 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2933 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2934 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2936 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2937 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2941 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2944 l_homo(k,ii)=.false.
2951 c Theta, dihedral and SC retraints
2953 if (waga_angle.gt.0.0d0) then
2954 c open (ientin,file=tpl_k_sigma_dih,status='old')
2955 c do irec=1,maxres-3 ! loop for reading sigma_dih
2956 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2957 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2958 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2959 c & sigma_dih(k,i+nnt-1)
2964 if (idomain(k,i).eq.0) then
2968 dih(k,i)=phiref(i) ! right?
2969 c read (ientin,*) sigma_dih(k,i) ! original variant
2970 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2971 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2972 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2973 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2974 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2976 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2977 & rescore(k,i-2)+rescore(k,i-3))/4.0
2978 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2979 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2980 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2981 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2982 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2983 c Instead of res sim other local measure of b/b str reliability possible
2984 if (sigma_dih(k,i).ne.0)
2985 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2986 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2991 if (waga_theta.gt.0.0d0) then
2992 c open (ientin,file=tpl_k_sigma_theta,status='old')
2993 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2994 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2995 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2996 c & sigma_theta(k,i+nnt-1)
3001 do i = nnt+2,nct ! right? without parallel.
3002 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3003 c do i=ithet_start,ithet_end ! with FG parallel.
3004 if (idomain(k,i).eq.0) then
3005 sigma_theta(k,i)=0.0
3008 thetatpl(k,i)=thetaref(i)
3009 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3010 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3011 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3012 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3013 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3014 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3015 & rescore(k,i-2))/3.0
3016 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3017 if (sigma_theta(k,i).ne.0)
3018 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3020 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3021 c rescore(k,i-2) ! right expression ?
3022 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3026 if (waga_d.gt.0.0d0) then
3027 c open (ientin,file=tpl_k_sigma_d,status='old')
3028 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3029 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3030 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3031 c & sigma_d(k,i+nnt-1)
3035 do i = nnt,nct ! right? without parallel.
3036 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3037 c do i=loc_start,loc_end ! with FG parallel.
3038 if (itype(i).eq.10) cycle
3039 if (idomain(k,i).eq.0 ) then
3046 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3047 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3048 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3049 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3050 sigma_d(k,i)=rescore(k,i) ! right expression ?
3051 if (sigma_d(k,i).ne.0)
3052 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3054 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3055 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3056 c read (ientin,*) sigma_d(k,i) ! 1st variant
3057 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3062 c remove distance restraints not used in any model from the list
3063 c shift data in all arrays
3065 if (waga_dist.ne.0.0d0) then
3071 if (ii_in_use(ii).eq.0.and.liiflag) then
3075 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3076 & .not.liiflag.and.ii.eq.lim_odl) then
3077 if (ii.eq.lim_odl) then
3078 iishift=ii-iistart+1
3083 do ki=iistart,lim_odl-iishift
3084 ires_homo(ki)=ires_homo(ki+iishift)
3085 jres_homo(ki)=jres_homo(ki+iishift)
3086 ii_in_use(ki)=ii_in_use(ki+iishift)
3087 do k=1,constr_homology
3088 odl(k,ki)=odl(k,ki+iishift)
3089 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3090 l_homo(k,ki)=l_homo(k,ki+iishift)
3094 lim_odl=lim_odl-iishift
3099 if (constr_homology.gt.0) call homology_partition
3100 if (constr_homology.gt.0) call init_int_table
3101 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3102 cd & "lim_xx=",lim_xx
3103 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3104 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3108 if (.not.lprn) return
3109 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3111 write (iout,*) "Distance restraints from templates"
3113 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3114 & ii,ires_homo(ii),jres_homo(ii),
3115 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3116 & ki=1,constr_homology)
3118 write (iout,*) "Dihedral angle restraints from templates"
3120 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3121 & (rad2deg*dih(ki,i),
3122 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3124 write (iout,*) "Virtual-bond angle restraints from templates"
3126 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3127 & (rad2deg*thetatpl(ki,i),
3128 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3130 write (iout,*) "SC restraints from templates"
3132 write(iout,'(i5,100(4f8.2,4x))') i,
3133 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3134 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3137 c -----------------------------------------------------------------
3140 c----------------------------------------------------------------------
3143 subroutine flush(iu)
3148 subroutine flush(iu)
3154 c------------------------------------------------------------------------------
3155 subroutine copy_to_tmp(source)
3156 include "DIMENSIONS"
3157 include "COMMON.IOUNITS"
3158 character*(*) source
3159 character* 256 tmpfile
3163 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3164 inquire(file=tmpfile,exist=ex)
3166 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3167 & " to temporary directory..."
3168 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3169 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3173 c------------------------------------------------------------------------------
3174 subroutine move_from_tmp(source)
3175 include "DIMENSIONS"
3176 include "COMMON.IOUNITS"
3177 character*(*) source
3180 write (*,*) "Moving ",source(:ilen(source)),
3181 & " from temporary directory to working directory"
3182 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3183 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3186 c------------------------------------------------------------------------------
3187 subroutine random_init(seed)
3189 C Initialize random number generator
3191 implicit real*8 (a-h,o-z)
3192 include 'DIMENSIONS'
3198 logical OKRandom, prng_restart
3200 integer iseed_array(4)
3202 include 'COMMON.IOUNITS'
3203 include 'COMMON.TIME1'
3204 include 'COMMON.THREAD'
3205 include 'COMMON.SBRIDGE'
3206 include 'COMMON.CONTROL'
3207 include 'COMMON.MCM'
3208 include 'COMMON.MAP'
3209 include 'COMMON.HEADER'
3210 csa include 'COMMON.CSA'
3211 include 'COMMON.CHAIN'
3212 include 'COMMON.MUCA'
3214 include 'COMMON.FFIELD'
3215 include 'COMMON.SETUP'
3216 iseed=-dint(dabs(seed))
3217 if (iseed.eq.0) then
3218 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3219 & 'Random seed undefined. The program will stop.'
3220 write (*,'(/80(1h*)/20x,a/80(1h*))')
3221 & 'Random seed undefined. The program will stop.'
3223 call mpi_finalize(mpi_comm_world,ierr)
3225 stop 'Bad random seed.'
3228 if (fg_rank.eq.0) then
3232 if(me.eq.king .or. .not. out1file)
3233 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3234 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3235 OKRandom = prng_restart(me,iseedi8)
3238 tmp=65536.0d0**(4-i)
3239 iseed_array(i) = dint(seed/tmp)
3240 seed=seed-iseed_array(i)*tmp
3242 if(me.eq.king .or. .not. out1file)
3243 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3244 & (iseed_array(i),i=1,4)
3245 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3246 & (iseed_array(i),i=1,4)
3247 OKRandom = prng_restart(me,iseed_array)
3250 r1=ran_number(0.0D0,1.0D0)
3251 if(me.eq.king .or. .not. out1file)
3252 & write (iout,*) 'ran_num',r1
3253 if (r1.lt.0.0d0) OKRandom=.false.
3255 if (.not.OKRandom) then
3256 write (iout,*) 'PRNG IS NOT WORKING!!!'
3257 print *,'PRNG IS NOT WORKING!!!'
3260 call mpi_abort(mpi_comm_world,error_msg,ierr)
3263 write (iout,*) 'too many processors for parallel prng'
3264 write (*,*) 'too many processors for parallel prng'
3272 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)