2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
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 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 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,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'IPRINT',iprint,0)
149 C SHIELD keyword sets if the shielding effect of side-chains is used
150 C 0 denots no shielding is used all peptide are equally despite the
151 C solvent accesible area
152 C 1 the newly introduced function
153 C 2 reseved for further possible developement
154 call readi(controlcard,'SHIELD',shield_mode,0)
155 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
156 write(iout,*) "shield_mode",shield_mode
158 call readi(controlcard,'MAXGEN',maxgen,10000)
159 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
160 call readi(controlcard,"KDIAG",kdiag,0)
161 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
162 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
163 & write (iout,*) "RESCALE_MODE",rescale_mode
164 split_ene=index(controlcard,'SPLIT_ENE').gt.0
165 if (index(controlcard,'REGULAR').gt.0.0D0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
170 if (index(controlcard,'CHECKGRAD').gt.0) then
172 if (index(controlcard,'CART').gt.0) then
174 elseif (index(controlcard,'CARINT').gt.0) then
179 elseif (index(controlcard,'THREAD').gt.0) then
181 call readi(controlcard,'THREAD',nthread,0)
182 if (nthread.gt.0) then
183 call reada(controlcard,'WEIDIS',weidis,0.1D0)
186 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
187 stop 'Error termination in Read_Control.'
189 else if (index(controlcard,'MCMA').gt.0) then
191 else if (index(controlcard,'MCEE').gt.0) then
193 else if (index(controlcard,'MULTCONF').gt.0) then
195 else if (index(controlcard,'MAP').gt.0) then
197 call readi(controlcard,'MAP',nmap,0)
198 else if (index(controlcard,'CSA').gt.0) then
200 crc else if (index(controlcard,'ZSCORE').gt.0) then
202 crc ZSCORE is rm from UNRES, modecalc=9 is available
205 cfcm else if (index(controlcard,'MCMF').gt.0) then
207 else if (index(controlcard,'SOFTREG').gt.0) then
209 else if (index(controlcard,'CHECK_BOND').gt.0) then
211 else if (index(controlcard,'TEST').gt.0) then
213 else if (index(controlcard,'MD').gt.0) then
215 else if (index(controlcard,'RE ').gt.0) then
219 lmuca=index(controlcard,'MUCA').gt.0
220 call readi(controlcard,'MUCADYN',mucadyn,0)
221 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
222 if (lmuca .and. (me.eq.king .or. .not.out1file ))
224 write (iout,*) 'MUCADYN=',mucadyn
225 write (iout,*) 'MUCASMOOTH=',muca_smooth
228 iscode=index(controlcard,'ONE_LETTER')
229 indphi=index(controlcard,'PHI')
230 indback=index(controlcard,'BACK')
231 iranconf=index(controlcard,'RAND_CONF')
232 i2ndstr=index(controlcard,'USE_SEC_PRED')
233 gradout=index(controlcard,'GRADOUT').gt.0
234 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
235 C DISTCHAINMAX become obsolete for periodic boundry condition
236 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
237 C Reading the dimensions of box in x,y,z coordinates
238 call reada(controlcard,'BOXX',boxxsize,100.0d0)
239 call reada(controlcard,'BOXY',boxysize,100.0d0)
240 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
241 c Cutoff range for interactions
242 call reada(controlcard,"R_CUT",r_cut,15.0d0)
243 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
244 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
245 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
246 if (lipthick.gt.0.0d0) then
247 bordliptop=(boxzsize+lipthick)/2.0
248 bordlipbot=bordliptop-lipthick
250 if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
251 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
252 buflipbot=bordlipbot+lipbufthick
253 bufliptop=bordliptop-lipbufthick
254 if ((lipbufthick*2.0d0).gt.lipthick)
255 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
257 write(iout,*) "bordliptop=",bordliptop
258 write(iout,*) "bordlipbot=",bordlipbot
259 write(iout,*) "bufliptop=",bufliptop
260 write(iout,*) "buflipbot=",buflipbot
261 write (iout,*) "SHIELD MODE",shield_mode
262 if (shield_mode.gt.0) then
264 C VSolvSphere the volume of solving sphere
266 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
267 C there will be no distinction between proline peptide group and normal peptide
268 C group in case of shielding parameters
269 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
270 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
271 write (iout,*) VSolvSphere,VSolvSphere_div
272 C long axis of side chain
274 long_r_sidechain(i)=vbldsc0(1,i)
275 short_r_sidechain(i)=sigma0(i)
279 if (me.eq.king .or. .not.out1file )
280 & write (iout,*) "DISTCHAINMAX",distchainmax
282 if(me.eq.king.or..not.out1file)
283 & write (iout,'(2a)') diagmeth(kdiag),
284 & ' routine used to diagonalize matrices.'
287 c--------------------------------------------------------------------------
288 subroutine read_REMDpar
292 implicit real*8 (a-h,o-z)
294 include 'COMMON.IOUNITS'
295 include 'COMMON.TIME1'
298 include 'COMMON.LANGEVIN'
300 include 'COMMON.LANGEVIN.lang0'
302 include 'COMMON.INTERACT'
303 include 'COMMON.NAMES'
305 include 'COMMON.REMD'
306 include 'COMMON.CONTROL'
307 include 'COMMON.SETUP'
309 character*320 controlcard
310 character*3200 controlcard1
311 integer iremd_m_total
313 if(me.eq.king.or..not.out1file)
314 & write (iout,*) "REMD setup"
316 call card_concat(controlcard)
317 call readi(controlcard,"NREP",nrep,3)
318 call readi(controlcard,"NSTEX",nstex,1000)
319 call reada(controlcard,"RETMIN",retmin,10.0d0)
320 call reada(controlcard,"RETMAX",retmax,1000.0d0)
321 mremdsync=(index(controlcard,'SYNC').gt.0)
322 call readi(controlcard,"NSYN",i_sync_step,100)
323 restart1file=(index(controlcard,'REST1FILE').gt.0)
324 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
325 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
326 if(max_cache_traj_use.gt.max_cache_traj)
327 & max_cache_traj_use=max_cache_traj
328 if(me.eq.king.or..not.out1file) then
329 cd if (traj1file) then
330 crc caching is in testing - NTWX is not ignored
331 cd write (iout,*) "NTWX value is ignored"
332 cd write (iout,*) " trajectory is stored to one file by master"
333 cd write (iout,*) " before exchange at NSTEX intervals"
335 write (iout,*) "NREP= ",nrep
336 write (iout,*) "NSTEX= ",nstex
337 write (iout,*) "SYNC= ",mremdsync
338 write (iout,*) "NSYN= ",i_sync_step
339 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
342 if (index(controlcard,'TLIST').gt.0) then
344 call card_concat(controlcard1)
345 read(controlcard1,*) (remd_t(i),i=1,nrep)
346 if(me.eq.king.or..not.out1file)
347 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
350 if (index(controlcard,'MLIST').gt.0) then
352 call card_concat(controlcard1)
353 read(controlcard1,*) (remd_m(i),i=1,nrep)
354 if(me.eq.king.or..not.out1file) then
355 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
358 iremd_m_total=iremd_m_total+remd_m(i)
360 write (iout,*) 'Total number of replicas ',iremd_m_total
363 if(me.eq.king.or..not.out1file)
364 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
367 c--------------------------------------------------------------------------
368 subroutine read_MDpar
372 implicit real*8 (a-h,o-z)
374 include 'COMMON.IOUNITS'
375 include 'COMMON.TIME1'
378 include 'COMMON.LANGEVIN'
380 include 'COMMON.LANGEVIN.lang0'
382 include 'COMMON.INTERACT'
383 include 'COMMON.NAMES'
385 include 'COMMON.SETUP'
386 include 'COMMON.CONTROL'
387 include 'COMMON.SPLITELE'
389 character*320 controlcard
391 call card_concat(controlcard)
392 call readi(controlcard,"NSTEP",n_timestep,1000000)
393 call readi(controlcard,"NTWE",ntwe,100)
394 call readi(controlcard,"NTWX",ntwx,1000)
395 call reada(controlcard,"DT",d_time,1.0d-1)
396 call reada(controlcard,"DVMAX",dvmax,2.0d1)
397 call reada(controlcard,"DAMAX",damax,1.0d1)
398 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
399 call readi(controlcard,"LANG",lang,0)
400 RESPA = index(controlcard,"RESPA") .gt. 0
401 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
402 ntime_split0=ntime_split
403 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
404 ntime_split0=ntime_split
405 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
406 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
407 rest = index(controlcard,"REST").gt.0
408 tbf = index(controlcard,"TBF").gt.0
409 usampl = index(controlcard,"USAMPL").gt.0
411 mdpdb = index(controlcard,"MDPDB").gt.0
412 call reada(controlcard,"T_BATH",t_bath,300.0d0)
413 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
414 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
415 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
416 if (count_reset_moment.eq.0) count_reset_moment=1000000000
417 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
418 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
419 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
420 if (count_reset_vel.eq.0) count_reset_vel=1000000000
421 large = index(controlcard,"LARGE").gt.0
422 print_compon = index(controlcard,"PRINT_COMPON").gt.0
423 rattle = index(controlcard,"RATTLE").gt.0
424 c if performing umbrella sampling, fragments constrained are read from the fragment file
430 if(me.eq.king.or..not.out1file) then
432 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
434 write (iout,'(a)') "The units are:"
435 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
436 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
437 & " acceleration: angstrom/(48.9 fs)**2"
438 write (iout,'(a)') "energy: kcal/mol, temperature: K"
440 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
441 write (iout,'(a60,f10.5,a)')
442 & "Initial time step of numerical integration:",d_time,
444 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
446 write (iout,'(2a,i4,a)')
447 & "A-MTS algorithm used; initial time step for fast-varying",
448 & " short-range forces split into",ntime_split," steps."
449 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
450 & r_cut," lambda",rlamb
452 write (iout,'(2a,f10.5)')
453 & "Maximum acceleration threshold to reduce the time step",
454 & "/increase split number:",damax
455 write (iout,'(2a,f10.5)')
456 & "Maximum predicted energy drift to reduce the timestep",
457 & "/increase split number:",edriftmax
458 write (iout,'(a60,f10.5)')
459 & "Maximum velocity threshold to reduce velocities:",dvmax
460 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
461 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
462 if (rattle) write (iout,'(a60)')
463 & "Rattle algorithm used to constrain the virtual bonds"
467 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
468 call reada(controlcard,"RWAT",rwat,1.4d0)
469 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
470 surfarea=index(controlcard,"SURFAREA").gt.0
471 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
472 if(me.eq.king.or..not.out1file)then
473 write (iout,'(/a,$)') "Langevin dynamics calculation"
476 & " with direct integration of Langevin equations"
477 else if (lang.eq.2) then
478 write (iout,'(a/)') " with TINKER stochasic MD integrator"
479 else if (lang.eq.3) then
480 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
481 else if (lang.eq.4) then
482 write (iout,'(a/)') " in overdamped mode"
484 write (iout,'(//a,i5)')
485 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
488 write (iout,'(a60,f10.5)') "Temperature:",t_bath
489 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
490 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
491 write (iout,'(a60,f10.5)')
492 & "Scaling factor of the friction forces:",scal_fric
493 if (surfarea) write (iout,'(2a,i10,a)')
494 & "Friction coefficients will be scaled by solvent-accessible",
495 & " surface area every",reset_fricmat," steps."
497 c Calculate friction coefficients and bounds of stochastic forces
498 eta=6*pi*cPoise*etawat
499 if(me.eq.king.or..not.out1file)
500 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
502 gamp=scal_fric*(pstok+rwat)*eta
503 stdfp=dsqrt(2*Rb*t_bath/d_time)
505 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
506 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
508 if(me.eq.king.or..not.out1file)then
509 write (iout,'(/2a/)')
510 & "Radii of site types and friction coefficients and std's of",
511 & " stochastic forces of fully exposed sites"
512 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
514 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
515 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
519 if(me.eq.king.or..not.out1file)then
520 write (iout,'(a)') "Berendsen bath calculation"
521 write (iout,'(a60,f10.5)') "Temperature:",t_bath
522 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
524 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
525 & count_reset_moment," steps"
527 & write (iout,'(a,i10,a)')
528 & "Velocities will be reset at random every",count_reset_vel,
532 if(me.eq.king.or..not.out1file)
533 & write (iout,'(a31)') "Microcanonical mode calculation"
535 if(me.eq.king.or..not.out1file)then
536 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
538 write(iout,*) "MD running with constraints."
539 write(iout,*) "Equilibration time ", eq_time, " mtus."
540 write(iout,*) "Constraining ", nfrag," fragments."
541 write(iout,*) "Length of each fragment, weight and q0:"
543 write (iout,*) "Set of restraints #",iset
545 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
546 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
548 write(iout,*) "constraints between ", npair, "fragments."
549 write(iout,*) "constraint pairs, weights and q0:"
551 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
552 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
554 write(iout,*) "angle constraints within ", nfrag_back,
555 & "backbone fragments."
556 write(iout,*) "fragment, weights:"
558 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
559 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
560 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
563 iset=mod(kolor,nset)+1
566 if(me.eq.king.or..not.out1file)
567 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
570 c------------------------------------------------------------------------------
573 C Read molecular data.
575 implicit real*8 (a-h,o-z)
581 include 'COMMON.IOUNITS'
584 include 'COMMON.INTERACT'
585 include 'COMMON.LOCAL'
586 include 'COMMON.NAMES'
587 include 'COMMON.CHAIN'
588 include 'COMMON.FFIELD'
589 include 'COMMON.SBRIDGE'
590 include 'COMMON.HEADER'
591 include 'COMMON.CONTROL'
592 include 'COMMON.DBASE'
593 include 'COMMON.THREAD'
594 include 'COMMON.CONTACTS'
595 include 'COMMON.TORCNSTR'
596 include 'COMMON.TIME1'
597 include 'COMMON.BOUNDS'
599 include 'COMMON.SETUP'
600 character*4 sequence(maxres)
602 double precision x(maxvar)
603 character*256 pdbfile
604 character*400 weightcard
605 character*80 weightcard_t,ucase
606 dimension itype_pdb(maxres)
607 common /pizda/ itype_pdb
608 logical seq_comp,fail
609 double precision energia(0:n_ene)
615 C Read weights of the subsequent energy terms.
616 call card_concat(weightcard)
617 call reada(weightcard,'WLONG',wlong,1.0D0)
618 call reada(weightcard,'WSC',wsc,wlong)
619 call reada(weightcard,'WSCP',wscp,wlong)
620 call reada(weightcard,'WELEC',welec,1.0D0)
621 call reada(weightcard,'WVDWPP',wvdwpp,welec)
622 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
623 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
624 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
625 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
626 call reada(weightcard,'WTURN3',wturn3,1.0D0)
627 call reada(weightcard,'WTURN4',wturn4,1.0D0)
628 call reada(weightcard,'WTURN6',wturn6,1.0D0)
629 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
630 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
631 call reada(weightcard,'WBOND',wbond,1.0D0)
632 call reada(weightcard,'WTOR',wtor,1.0D0)
633 call reada(weightcard,'WTORD',wtor_d,1.0D0)
634 call reada(weightcard,'WANG',wang,1.0D0)
635 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
636 call reada(weightcard,'SCAL14',scal14,0.4D0)
637 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
638 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
639 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
640 call reada(weightcard,'TEMP0',temp0,300.0d0)
641 call reada(weightcard,'WLT',wliptran,0.0D0)
642 if (index(weightcard,'SOFT').gt.0) ipot=6
643 C 12/1/95 Added weight for the multi-body term WCORR
644 call reada(weightcard,'WCORRH',wcorr,1.0D0)
645 if (wcorr4.gt.0.0d0) wcorr=wcorr4
665 if(me.eq.king.or..not.out1file)
666 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
667 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
669 10 format (/'Energy-term weights (unscaled):'//
670 & 'WSCC= ',f10.6,' (SC-SC)'/
671 & 'WSCP= ',f10.6,' (SC-p)'/
672 & 'WELEC= ',f10.6,' (p-p electr)'/
673 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
674 & 'WBOND= ',f10.6,' (stretching)'/
675 & 'WANG= ',f10.6,' (bending)'/
676 & 'WSCLOC= ',f10.6,' (SC local)'/
677 & 'WTOR= ',f10.6,' (torsional)'/
678 & 'WTORD= ',f10.6,' (double torsional)'/
679 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
680 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
681 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
682 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
683 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
684 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
685 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
686 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
687 & 'WTURN6= ',f10.6,' (turns, 6th order)')
688 if(me.eq.king.or..not.out1file)then
689 if (wcorr4.gt.0.0d0) then
690 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
691 & 'between contact pairs of peptide groups'
692 write (iout,'(2(a,f5.3/))')
693 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
694 & 'Range of quenching the correlation terms:',2*delt_corr
695 else if (wcorr.gt.0.0d0) then
696 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
697 & 'between contact pairs of peptide groups'
699 write (iout,'(a,f8.3)')
700 & 'Scaling factor of 1,4 SC-p interactions:',scal14
701 write (iout,'(a,f8.3)')
702 & 'General scaling factor of SC-p interactions:',scalscp
704 r0_corr=cutoff_corr-delt_corr
706 aad(i,1)=scalscp*aad(i,1)
707 aad(i,2)=scalscp*aad(i,2)
708 bad(i,1)=scalscp*bad(i,1)
709 bad(i,2)=scalscp*bad(i,2)
711 call rescale_weights(t_bath)
712 if(me.eq.king.or..not.out1file)
713 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
714 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
716 22 format (/'Energy-term weights (scaled):'//
717 & 'WSCC= ',f10.6,' (SC-SC)'/
718 & 'WSCP= ',f10.6,' (SC-p)'/
719 & 'WELEC= ',f10.6,' (p-p electr)'/
720 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
721 & 'WBOND= ',f10.6,' (stretching)'/
722 & 'WANG= ',f10.6,' (bending)'/
723 & 'WSCLOC= ',f10.6,' (SC local)'/
724 & 'WTOR= ',f10.6,' (torsional)'/
725 & 'WTORD= ',f10.6,' (double torsional)'/
726 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
727 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
728 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
729 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
730 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
731 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
732 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
733 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
734 & 'WTURN6= ',f10.6,' (turns, 6th order)')
735 if(me.eq.king.or..not.out1file)
736 & write (iout,*) "Reference temperature for weights calculation:",
738 call reada(weightcard,"D0CM",d0cm,3.78d0)
739 call reada(weightcard,"AKCM",akcm,15.1d0)
740 call reada(weightcard,"AKTH",akth,11.0d0)
741 call reada(weightcard,"AKCT",akct,12.0d0)
742 call reada(weightcard,"V1SS",v1ss,-1.08d0)
743 call reada(weightcard,"V2SS",v2ss,7.61d0)
744 call reada(weightcard,"V3SS",v3ss,13.7d0)
745 call reada(weightcard,"EBR",ebr,-5.50D0)
746 call reada(weightcard,"ATRISS",atriss,0.301D0)
747 call reada(weightcard,"BTRISS",btriss,0.021D0)
748 call reada(weightcard,"CTRISS",ctriss,1.001D0)
749 call reada(weightcard,"DTRISS",dtriss,1.001D0)
750 write (iout,*) "ATRISS=", atriss
751 write (iout,*) "BTRISS=", btriss
752 write (iout,*) "CTRISS=", ctriss
753 write (iout,*) "DTRISS=", dtriss
754 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
756 dyn_ss_mask(i)=.false.
760 dyn_ssbond_ij(i,j)=1.0d300
763 call reada(weightcard,"HT",Ht,0.0D0)
765 ss_depth=ebr/wsc-0.25*eps(1,1)
766 Ht=Ht/wsc-0.25*eps(1,1)
767 akcm=akcm*wstrain/wsc
768 akth=akth*wstrain/wsc
769 akct=akct*wstrain/wsc
770 v1ss=v1ss*wstrain/wsc
771 v2ss=v2ss*wstrain/wsc
772 v3ss=v3ss*wstrain/wsc
774 if (wstrain.ne.0.0) then
775 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
781 if(me.eq.king.or..not.out1file) then
782 write (iout,*) "Parameters of the SS-bond potential:"
783 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
785 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
786 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
787 write (iout,*)" HT",Ht
788 print *,'indpdb=',indpdb,' pdbref=',pdbref
790 if (indpdb.gt.0 .or. pdbref) then
791 read(inp,'(a)') pdbfile
792 if(me.eq.king.or..not.out1file)
793 & write (iout,'(2a)') 'PDB data will be read from file ',
794 & pdbfile(:ilen(pdbfile))
795 open(ipdbin,file=pdbfile,status='old',err=33)
797 33 write (iout,'(a)') 'Error opening PDB file.'
800 c write (iout,*) 'Begin reading pdb data'
803 c write (iout,*) 'Finished reading pdb data'
805 if(me.eq.king.or..not.out1file)
806 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
807 & ' nstart_sup=',nstart_sup
809 itype_pdb(i)=itype(i)
813 nct=nstart_sup+nsup-1
814 call contact(.false.,ncont_ref,icont_ref,co)
817 if(me.eq.king.or..not.out1file)
818 & write(iout,*)'Adding sidechains'
822 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
825 do while (fail.and.nsi.le.maxsi)
826 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
829 if(fail) write(iout,*)'Adding sidechain failed for res ',
830 & i,' after ',nsi,' trials'
835 if (indpdb.eq.0) then
836 C Read sequence if not taken from the pdb file.
838 c print *,'nres=',nres
839 if (iscode.gt.0) then
840 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
842 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
844 C Convert sequence to numeric code
846 itype(i)=rescode(i,sequence(i),iscode)
848 C Assign initial virtual bond lengths
854 vbld(i+nres)=dsc(iabs(itype(i)))
855 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
856 c write (iout,*) "i",i," itype",itype(i),
857 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
861 c print '(20i4)',(itype(i),i=1,nres)
864 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
866 if (itype(i).eq.ntyp1) then
870 else if (iabs(itype(i+1)).ne.20) then
872 else if (iabs(itype(i)).ne.20) then
879 if(me.eq.king.or..not.out1file)then
880 write (iout,*) "ITEL"
882 write (iout,*) i,itype(i),itel(i)
884 print *,'Call Read_Bridge.'
887 C 8/13/98 Set limits to generating the dihedral angles
892 read (inp,*) ndih_constr
893 if (ndih_constr.gt.0) then
895 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
897 if(me.eq.king.or..not.out1file)then
899 & 'There are',ndih_constr,' constraints on phi angles.'
901 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
906 phi0(i)=deg2rad*phi0(i)
907 drange(i)=deg2rad*drange(i)
909 C if(me.eq.king.or..not.out1file)
910 C & write (iout,*) 'FTORS',ftors
913 phibound(1,ii) = phi0(i)-drange(i)
914 phibound(2,ii) = phi0(i)+drange(i)
917 C first setting the theta boundaries to 0 to pi
918 C this mean that there is no energy penalty for any angle occuring this can be applied
919 C for generate random conformation but is not implemented in this way
924 C begin reading theta constrains this is quartic constrains allowing to
925 C have smooth second derivative
926 if (with_theta_constr) then
927 C with_theta_constr is keyword allowing for occurance of theta constrains
928 read (inp,*) ntheta_constr
929 C ntheta_constr is the number of theta constrains
930 if (ntheta_constr.gt.0) then
932 read (inp,*) (itheta_constr(i),theta_constr0(i),
933 & theta_drange(i),for_thet_constr(i),
935 C the above code reads from 1 to ntheta_constr
936 C itheta_constr(i) residue i for which is theta_constr
937 C theta_constr0 the global minimum value
938 C theta_drange is range for which there is no energy penalty
939 C for_thet_constr is the force constant for quartic energy penalty
941 if(me.eq.king.or..not.out1file)then
943 & 'There are',ntheta_constr,' constraints on phi angles.'
945 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
951 theta_constr0(i)=deg2rad*theta_constr0(i)
952 theta_drange(i)=deg2rad*theta_drange(i)
954 C if(me.eq.king.or..not.out1file)
955 C & write (iout,*) 'FTORS',ftors
956 C do i=1,ntheta_constr
957 C ii = itheta_constr(i)
958 C thetabound(1,ii) = phi0(i)-drange(i)
959 C thetabound(2,ii) = phi0(i)+drange(i)
961 endif ! ntheta_constr.gt.0
962 endif! with_theta_constr
964 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
965 C write (iout,*) "with_dihed_constr ",with_dihed_constr
970 write (iout,'(a)') 'Boundaries in phi angle sampling:'
972 write (iout,'(a3,i5,2f10.1)')
973 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
979 cd print *,'NNT=',NNT,' NCT=',NCT
980 if (itype(1).eq.ntyp1) nnt=2
981 if (itype(nres).eq.ntyp1) nct=nct-1
983 if(me.eq.king.or..not.out1file)
984 & write (iout,'(a,i3)') 'nsup=',nsup
986 if (nsup.le.(nct-nnt+1)) then
987 do i=0,nct-nnt+1-nsup
988 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
994 & 'Error - sequences to be superposed do not match.'
997 do i=0,nsup-(nct-nnt+1)
998 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1000 nstart_sup=nstart_sup+i
1006 & 'Error - sequences to be superposed do not match.'
1009 if (nsup.eq.0) nsup=nct-nnt
1010 if (nstart_sup.eq.0) nstart_sup=nnt
1011 if (nstart_seq.eq.0) nstart_seq=nnt
1012 if(me.eq.king.or..not.out1file)
1013 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1014 & ' nstart_seq=',nstart_seq
1016 c--- Zscore rms -------
1017 if (nz_start.eq.0) nz_start=nnt
1018 if (nz_end.eq.0 .and. nsup.gt.0) then
1020 else if (nz_end.eq.0) then
1023 if(me.eq.king.or..not.out1file)then
1024 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1025 write (iout,*) 'IZ_SC=',iz_sc
1027 c----------------------
1030 if (.not.pdbref) then
1031 call read_angles(inp,*38)
1033 38 write (iout,'(a)') 'Error reading reference structure.'
1035 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1036 stop 'Error reading reference structure'
1040 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1047 cref(j,i,kkk)=c(j,i)
1050 call contact(.true.,ncont_ref,icont_ref,co)
1054 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1056 if (constr_dist.gt.0) call read_dist_constr
1057 write (iout,*) "After read_dist_constr nhpb",nhpb
1058 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1060 if(me.eq.king.or..not.out1file)
1061 & write (iout,*) 'Contact order:',co
1063 if(me.eq.king.or..not.out1file)
1064 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1067 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1069 if(me.eq.king.or..not.out1file)
1070 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1071 & icont_ref(1,i),' ',
1072 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1076 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1077 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1078 & modecalc.ne.10) then
1079 C If input structure hasn't been supplied from the PDB file read or generate
1081 if (iranconf.eq.0 .and. .not. extconf) then
1082 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1083 & write (iout,'(a)') 'Initial geometry will be read in.'
1085 read(inp,'(8f10.5)',end=36,err=36)
1086 & ((c(l,k),l=1,3),k=1,nres),
1087 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1088 write (iout,*) "Exit READ_CART"
1089 write (iout,'(8f10.5)')
1090 & ((c(l,k),l=1,3),k=1,nres),
1091 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1092 call int_from_cart1(.true.)
1093 write (iout,*) "Finish INT_TO_CART"
1096 dc(j,i)=c(j,i+1)-c(j,i)
1097 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1101 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1103 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1104 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1110 call read_angles(inp,*36)
1113 36 write (iout,'(a)') 'Error reading angle file.'
1115 call mpi_finalize( MPI_COMM_WORLD,IERR )
1117 stop 'Error reading angle file.'
1119 else if (extconf) then
1120 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1121 & write (iout,'(a)') 'Extended chain initial geometry.'
1123 theta(i)=90d0*deg2rad
1126 phi(i)=180d0*deg2rad
1129 alph(i)=110d0*deg2rad
1132 omeg(i)=-120d0*deg2rad
1133 if (itype(i).le.0) omeg(i)=-omeg(i)
1136 if(me.eq.king.or..not.out1file)
1137 & write (iout,'(a)') 'Random-generated initial geometry.'
1141 if (me.eq.king .or. fg_rank.eq.0 .and. (
1142 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1146 call gen_rand_conf(itmp,*30)
1148 30 write (iout,*) 'Failed to generate random conformation',
1149 & ', itrial=',itrial
1150 write (*,*) 'Processor:',me,
1151 & ' Failed to generate random conformation',
1161 write (iout,'(a,i3,a)') 'Processor:',me,
1162 & ' error in generating random conformation.'
1163 write (*,'(a,i3,a)') 'Processor:',me,
1164 & ' error in generating random conformation.'
1167 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1172 & ' error in generating random conformation.'
1177 elseif (modecalc.eq.4) then
1178 read (inp,'(a)') intinname
1179 open (intin,file=intinname,status='old',err=333)
1180 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1181 & write (iout,'(a)') 'intinname',intinname
1182 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1184 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1186 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1188 stop 'Error opening angle file.'
1192 C Generate distance constraints, if the PDB structure is to be regularized.
1193 if (nthread.gt.0) then
1194 call read_threadbase
1197 if (me.eq.king .or. .not. out1file)
1199 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1200 write (iout,'(/a,i3,a)')
1201 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1202 write (iout,'(20i4)') (iss(i),i=1,ns)
1204 write(iout,*)"Running with dynamic disulfide-bond formation"
1206 write (iout,'(/a/)') 'Pre-formed links are:'
1212 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1213 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1219 if (ns.gt.0.and.dyn_ss) then
1223 forcon(i-nss)=forcon(i)
1230 dyn_ss_mask(iss(i))=.true.
1233 if (i2ndstr.gt.0) call secstrp2dihc
1234 c call geom_to_var(nvar,x)
1235 c call etotal(energia(0))
1236 c call enerprint(energia(0))
1237 c call briefout(0,etot)
1239 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1240 cd write (iout,'(a)') 'Variable list:'
1241 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1243 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1244 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1245 & 'Processor',myrank,': end reading molecular data.'
1250 c--------------------------------------------------------------------------
1251 logical function seq_comp(itypea,itypeb,length)
1253 integer length,itypea(length),itypeb(length)
1256 if (itypea(i).ne.itypeb(i)) then
1264 c-----------------------------------------------------------------------------
1265 subroutine read_bridge
1266 C Read information about disulfide bridges.
1267 implicit real*8 (a-h,o-z)
1268 include 'DIMENSIONS'
1272 include 'COMMON.IOUNITS'
1273 include 'COMMON.GEO'
1274 include 'COMMON.VAR'
1275 include 'COMMON.INTERACT'
1276 include 'COMMON.LOCAL'
1277 include 'COMMON.NAMES'
1278 include 'COMMON.CHAIN'
1279 include 'COMMON.FFIELD'
1280 include 'COMMON.SBRIDGE'
1281 include 'COMMON.HEADER'
1282 include 'COMMON.CONTROL'
1283 include 'COMMON.DBASE'
1284 include 'COMMON.THREAD'
1285 include 'COMMON.TIME1'
1286 include 'COMMON.SETUP'
1287 C Read bridging residues.
1288 read (inp,*) ns,(iss(i),i=1,ns)
1290 if(me.eq.king.or..not.out1file)
1291 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1292 C Check whether the specified bridging residues are cystines.
1294 if (itype(iss(i)).ne.1) then
1295 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1296 & 'Do you REALLY think that the residue ',
1297 & restyp(itype(iss(i))),i,
1298 & ' can form a disulfide bridge?!!!'
1299 write (*,'(2a,i3,a)')
1300 & 'Do you REALLY think that the residue ',
1301 & restyp(itype(iss(i))),i,
1302 & ' can form a disulfide bridge?!!!'
1304 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1309 C Read preformed bridges.
1311 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1313 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1316 C Check if the residues involved in bridges are in the specified list of
1317 C bridging residues.
1320 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1321 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1322 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1323 & ' contains residues present in other pairs.'
1324 write (*,'(a,i3,a)') 'Disulfide pair',i,
1325 & ' contains residues present in other pairs.'
1327 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1333 if (ihpb(i).eq.iss(j)) goto 10
1335 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1338 if (jhpb(i).eq.iss(j)) goto 20
1340 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1346 ihpb(i)=ihpb(i)+nres
1347 jhpb(i)=jhpb(i)+nres
1353 c----------------------------------------------------------------------------
1354 subroutine read_x(kanal,*)
1355 implicit real*8 (a-h,o-z)
1356 include 'DIMENSIONS'
1357 include 'COMMON.GEO'
1358 include 'COMMON.VAR'
1359 include 'COMMON.CHAIN'
1360 include 'COMMON.IOUNITS'
1361 include 'COMMON.CONTROL'
1362 include 'COMMON.LOCAL'
1363 include 'COMMON.INTERACT'
1364 c Read coordinates from input
1366 read(kanal,'(8f10.5)',end=10,err=10)
1367 & ((c(l,k),l=1,3),k=1,nres),
1368 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1371 c(j,2*nres)=c(j,nres)
1373 call int_from_cart1(.false.)
1376 dc(j,i)=c(j,i+1)-c(j,i)
1377 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1381 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1383 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1384 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1392 c----------------------------------------------------------------------------
1393 subroutine read_threadbase
1394 implicit real*8 (a-h,o-z)
1395 include 'DIMENSIONS'
1396 include 'COMMON.IOUNITS'
1397 include 'COMMON.GEO'
1398 include 'COMMON.VAR'
1399 include 'COMMON.INTERACT'
1400 include 'COMMON.LOCAL'
1401 include 'COMMON.NAMES'
1402 include 'COMMON.CHAIN'
1403 include 'COMMON.FFIELD'
1404 include 'COMMON.SBRIDGE'
1405 include 'COMMON.HEADER'
1406 include 'COMMON.CONTROL'
1407 include 'COMMON.DBASE'
1408 include 'COMMON.THREAD'
1409 include 'COMMON.TIME1'
1410 C Read pattern database for threading.
1411 read (icbase,*) nseq
1413 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1414 & nres_base(2,i),nres_base(3,i)
1415 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1417 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1418 c & nres_base(2,i),nres_base(3,i)
1419 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1423 if (weidis.eq.0.0D0) weidis=0.1D0
1432 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1433 write (iout,'(a,i5)') 'nexcl: ',nexcl
1434 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1437 c------------------------------------------------------------------------------
1438 subroutine setup_var
1439 implicit real*8 (a-h,o-z)
1440 include 'DIMENSIONS'
1441 include 'COMMON.IOUNITS'
1442 include 'COMMON.GEO'
1443 include 'COMMON.VAR'
1444 include 'COMMON.INTERACT'
1445 include 'COMMON.LOCAL'
1446 include 'COMMON.NAMES'
1447 include 'COMMON.CHAIN'
1448 include 'COMMON.FFIELD'
1449 include 'COMMON.SBRIDGE'
1450 include 'COMMON.HEADER'
1451 include 'COMMON.CONTROL'
1452 include 'COMMON.DBASE'
1453 include 'COMMON.THREAD'
1454 include 'COMMON.TIME1'
1455 C Set up variable list.
1461 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1463 ialph(i,1)=nvar+nside
1467 if (indphi.gt.0) then
1469 else if (indback.gt.0) then
1474 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1477 c----------------------------------------------------------------------------
1478 subroutine gen_dist_constr
1479 C Generate CA distance constraints.
1480 implicit real*8 (a-h,o-z)
1481 include 'DIMENSIONS'
1482 include 'COMMON.IOUNITS'
1483 include 'COMMON.GEO'
1484 include 'COMMON.VAR'
1485 include 'COMMON.INTERACT'
1486 include 'COMMON.LOCAL'
1487 include 'COMMON.NAMES'
1488 include 'COMMON.CHAIN'
1489 include 'COMMON.FFIELD'
1490 include 'COMMON.SBRIDGE'
1491 include 'COMMON.HEADER'
1492 include 'COMMON.CONTROL'
1493 include 'COMMON.DBASE'
1494 include 'COMMON.THREAD'
1495 include 'COMMON.TIME1'
1496 dimension itype_pdb(maxres)
1497 common /pizda/ itype_pdb
1499 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1500 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1501 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1503 do i=nstart_sup,nstart_sup+nsup-1
1504 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1505 cd & ' seq_pdb', restyp(itype_pdb(i))
1506 do j=i+2,nstart_sup+nsup-1
1508 ihpb(nhpb)=i+nstart_seq-nstart_sup
1509 jhpb(nhpb)=j+nstart_seq-nstart_sup
1511 dhpb(nhpb)=dist(i,j)
1514 cd write (iout,'(a)') 'Distance constraints:'
1519 cd if (ii.gt.nres) then
1524 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1525 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1526 cd & dhpb(i),forcon(i)
1530 c----------------------------------------------------------------------------
1532 implicit real*8 (a-h,o-z)
1533 include 'DIMENSIONS'
1534 include 'COMMON.MAP'
1535 include 'COMMON.IOUNITS'
1536 character*3 angid(4) /'THE','PHI','ALP','OME'/
1537 character*80 mapcard,ucase
1539 read (inp,'(a)') mapcard
1540 mapcard=ucase(mapcard)
1541 if (index(mapcard,'PHI').gt.0) then
1543 else if (index(mapcard,'THE').gt.0) then
1545 else if (index(mapcard,'ALP').gt.0) then
1547 else if (index(mapcard,'OME').gt.0) then
1550 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1551 stop 'Error - illegal variable spec in MAP card.'
1553 call readi (mapcard,'RES1',res1(imap),0)
1554 call readi (mapcard,'RES2',res2(imap),0)
1555 if (res1(imap).eq.0) then
1556 res1(imap)=res2(imap)
1557 else if (res2(imap).eq.0) then
1558 res2(imap)=res1(imap)
1560 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1562 & 'Error - illegal definition of variable group in MAP.'
1563 stop 'Error - illegal definition of variable group in MAP.'
1565 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1566 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1567 call readi(mapcard,'NSTEP',nstep(imap),0)
1568 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1570 & 'Illegal boundary and/or step size specification in MAP.'
1571 stop 'Illegal boundary and/or step size specification in MAP.'
1576 c----------------------------------------------------------------------------
1578 implicit real*8 (a-h,o-z)
1579 include 'DIMENSIONS'
1580 include 'COMMON.IOUNITS'
1581 include 'COMMON.GEO'
1582 include 'COMMON.CSA'
1583 include 'COMMON.BANK'
1584 include 'COMMON.CONTROL'
1586 character*620 mcmcard
1587 call card_concat(mcmcard)
1589 call readi(mcmcard,'NCONF',nconf,50)
1590 call readi(mcmcard,'NADD',nadd,0)
1591 call readi(mcmcard,'JSTART',jstart,1)
1592 call readi(mcmcard,'JEND',jend,1)
1593 call readi(mcmcard,'NSTMAX',nstmax,500000)
1594 call readi(mcmcard,'N0',n0,1)
1595 call readi(mcmcard,'N1',n1,6)
1596 call readi(mcmcard,'N2',n2,4)
1597 call readi(mcmcard,'N3',n3,0)
1598 call readi(mcmcard,'N4',n4,0)
1599 call readi(mcmcard,'N5',n5,0)
1600 call readi(mcmcard,'N6',n6,10)
1601 call readi(mcmcard,'N7',n7,0)
1602 call readi(mcmcard,'N8',n8,0)
1603 call readi(mcmcard,'N9',n9,0)
1604 call readi(mcmcard,'N14',n14,0)
1605 call readi(mcmcard,'N15',n15,0)
1606 call readi(mcmcard,'N16',n16,0)
1607 call readi(mcmcard,'N17',n17,0)
1608 call readi(mcmcard,'N18',n18,0)
1610 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1612 call readi(mcmcard,'NDIFF',ndiff,2)
1613 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1614 call readi(mcmcard,'IS1',is1,1)
1615 call readi(mcmcard,'IS2',is2,8)
1616 call readi(mcmcard,'NRAN0',nran0,4)
1617 call readi(mcmcard,'NRAN1',nran1,2)
1618 call readi(mcmcard,'IRR',irr,1)
1619 call readi(mcmcard,'NSEED',nseed,20)
1620 call readi(mcmcard,'NTOTAL',ntotal,10000)
1621 call reada(mcmcard,'CUT1',cut1,2.0d0)
1622 call reada(mcmcard,'CUT2',cut2,5.0d0)
1623 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1624 call readi(mcmcard,'ICMAX',icmax,3)
1625 call readi(mcmcard,'IRESTART',irestart,0)
1626 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1629 call reada(mcmcard,'DELE',dele,20.0d0)
1630 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1631 call readi(mcmcard,'IREF',iref,0)
1632 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1633 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1634 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1635 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1636 write (iout,*) "NCONF_IN",nconf_in
1639 c----------------------------------------------------------------------------
1640 cfmc subroutine mcmfread
1641 cfmc implicit real*8 (a-h,o-z)
1642 cfmc include 'DIMENSIONS'
1643 cfmc include 'COMMON.MCMF'
1644 cfmc include 'COMMON.IOUNITS'
1645 cfmc include 'COMMON.GEO'
1646 cfmc character*80 ucase
1647 cfmc character*620 mcmcard
1648 cfmc call card_concat(mcmcard)
1650 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1651 cfmc write(iout,*)'MAXRANT=',maxrant
1652 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1653 cfmc write(iout,*)'MAXFAM=',maxfam
1654 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1655 cfmc write(iout,*)'NNET1=',nnet1
1656 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1657 cfmc write(iout,*)'NNET2=',nnet2
1658 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1659 cfmc write(iout,*)'NNET3=',nnet3
1660 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1661 cfmc write(iout,*)'ILASTT=',ilastt
1662 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1663 cfmc write(iout,*)'MAXSTR=',maxstr
1664 cfmc maxstr_f=maxstr/maxfam
1665 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1666 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1667 cfmc write(iout,*)'NMCMF=',nmcmf
1668 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1669 cfmc write(iout,*)'IFOCUS=',ifocus
1670 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1671 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1672 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1673 cfmc write(iout,*)'INTPRT=',intprt
1674 cfmc call readi(mcmcard,'IPRT',iprt,100)
1675 cfmc write(iout,*)'IPRT=',iprt
1676 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1677 cfmc write(iout,*)'IMAXTR=',imaxtr
1678 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1679 cfmc write(iout,*)'MAXEVEN=',maxeven
1680 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1681 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1682 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1683 cfmc write(iout,*)'INIMIN=',inimin
1684 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1685 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1686 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1687 cfmc write(iout,*)'NTHREAD=',nthread
1688 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1689 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1690 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1691 cfmc write(iout,*)'MAXPERT=',maxpert
1692 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1693 cfmc write(iout,*)'IRMSD=',irmsd
1694 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1695 cfmc write(iout,*)'DENEMIN=',denemin
1696 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1697 cfmc write(iout,*)'RCUT1S=',rcut1s
1698 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1699 cfmc write(iout,*)'RCUT1E=',rcut1e
1700 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1701 cfmc write(iout,*)'RCUT2S=',rcut2s
1702 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1703 cfmc write(iout,*)'RCUT2E=',rcut2e
1704 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1705 cfmc write(iout,*)'DPERT1=',d_pert1
1706 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1707 cfmc write(iout,*)'DPERT1A=',d_pert1a
1708 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1709 cfmc write(iout,*)'DPERT2=',d_pert2
1710 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1711 cfmc write(iout,*)'DPERT2A=',d_pert2a
1712 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1713 cfmc write(iout,*)'DPERT2B=',d_pert2b
1714 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1715 cfmc write(iout,*)'DPERT2C=',d_pert2c
1716 cfmc d_pert1=deg2rad*d_pert1
1717 cfmc d_pert1a=deg2rad*d_pert1a
1718 cfmc d_pert2=deg2rad*d_pert2
1719 cfmc d_pert2a=deg2rad*d_pert2a
1720 cfmc d_pert2b=deg2rad*d_pert2b
1721 cfmc d_pert2c=deg2rad*d_pert2c
1722 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1723 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1724 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1725 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1726 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1727 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1728 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1729 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1730 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1731 cfmc write(iout,*)'RCUTINI=',rcutini
1732 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1733 cfmc write(iout,*)'GRAT=',grat
1734 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1735 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1739 c----------------------------------------------------------------------------
1741 implicit real*8 (a-h,o-z)
1742 include 'DIMENSIONS'
1743 include 'COMMON.MCM'
1744 include 'COMMON.MCE'
1745 include 'COMMON.IOUNITS'
1747 character*320 mcmcard
1748 call card_concat(mcmcard)
1749 call readi(mcmcard,'MAXACC',maxacc,100)
1750 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1751 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1752 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1753 call readi(mcmcard,'MAXREPM',maxrepm,200)
1754 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1755 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1756 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1757 call reada(mcmcard,'E_UP',e_up,5.0D0)
1758 call reada(mcmcard,'DELTE',delte,0.1D0)
1759 call readi(mcmcard,'NSWEEP',nsweep,5)
1760 call readi(mcmcard,'NSTEPH',nsteph,0)
1761 call readi(mcmcard,'NSTEPC',nstepc,0)
1762 call reada(mcmcard,'TMIN',tmin,298.0D0)
1763 call reada(mcmcard,'TMAX',tmax,298.0D0)
1764 call readi(mcmcard,'NWINDOW',nwindow,0)
1765 call readi(mcmcard,'PRINT_MC',print_mc,0)
1766 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1767 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1768 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1769 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1770 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1771 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1772 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1773 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1774 if (nwindow.gt.0) then
1775 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1777 winlen(i)=winend(i)-winstart(i)+1
1780 if (tmax.lt.tmin) tmax=tmin
1781 if (tmax.eq.tmin) then
1785 if (nstepc.gt.0 .and. nsteph.gt.0) then
1786 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1787 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1789 C Probabilities of different move types
1790 sumpro_type(0)=0.0D0
1791 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1792 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1793 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1794 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1795 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1796 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1797 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1799 print *,'i',i,' sumprotype',sumpro_type(i)
1800 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1801 print *,'i',i,' sumprotype',sumpro_type(i)
1805 c----------------------------------------------------------------------------
1806 subroutine read_minim
1807 implicit real*8 (a-h,o-z)
1808 include 'DIMENSIONS'
1809 include 'COMMON.MINIM'
1810 include 'COMMON.IOUNITS'
1812 character*320 minimcard
1813 call card_concat(minimcard)
1814 call readi(minimcard,'MAXMIN',maxmin,2000)
1815 call readi(minimcard,'MAXFUN',maxfun,5000)
1816 call readi(minimcard,'MINMIN',minmin,maxmin)
1817 call readi(minimcard,'MINFUN',minfun,maxmin)
1818 call reada(minimcard,'TOLF',tolf,1.0D-2)
1819 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1820 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1821 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1822 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1823 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1824 & 'Options in energy minimization:'
1825 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1826 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1827 & 'MinMin:',MinMin,' MinFun:',MinFun,
1828 & ' TolF:',TolF,' RTolF:',RTolF
1831 c----------------------------------------------------------------------------
1832 subroutine read_angles(kanal,*)
1833 implicit real*8 (a-h,o-z)
1834 include 'DIMENSIONS'
1835 include 'COMMON.GEO'
1836 include 'COMMON.VAR'
1837 include 'COMMON.CHAIN'
1838 include 'COMMON.IOUNITS'
1839 include 'COMMON.CONTROL'
1840 c Read angles from input
1842 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1843 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1844 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1845 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1848 c 9/7/01 avoid 180 deg valence angle
1849 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1851 theta(i)=deg2rad*theta(i)
1852 phi(i)=deg2rad*phi(i)
1853 alph(i)=deg2rad*alph(i)
1854 omeg(i)=deg2rad*omeg(i)
1859 c----------------------------------------------------------------------------
1860 subroutine reada(rekord,lancuch,wartosc,default)
1862 character*(*) rekord,lancuch
1863 double precision wartosc,default
1866 iread=index(rekord,lancuch)
1867 if (iread.eq.0) then
1871 iread=iread+ilen(lancuch)+1
1872 read (rekord(iread:),*,err=10,end=10) wartosc
1877 c----------------------------------------------------------------------------
1878 subroutine readi(rekord,lancuch,wartosc,default)
1880 character*(*) rekord,lancuch
1881 integer wartosc,default
1884 iread=index(rekord,lancuch)
1885 if (iread.eq.0) then
1889 iread=iread+ilen(lancuch)+1
1890 read (rekord(iread:),*,err=10,end=10) wartosc
1895 c----------------------------------------------------------------------------
1896 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1899 integer tablica(dim),default
1900 character*(*) rekord,lancuch
1907 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1908 if (iread.eq.0) return
1909 iread=iread+ilen(lancuch)+1
1910 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1913 c----------------------------------------------------------------------------
1914 subroutine multreada(rekord,lancuch,tablica,dim,default)
1917 double precision tablica(dim),default
1918 character*(*) rekord,lancuch
1925 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1926 if (iread.eq.0) return
1927 iread=iread+ilen(lancuch)+1
1928 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1931 c----------------------------------------------------------------------------
1932 subroutine openunits
1933 implicit real*8 (a-h,o-z)
1934 include 'DIMENSIONS'
1937 character*16 form,nodename
1940 include 'COMMON.SETUP'
1941 include 'COMMON.IOUNITS'
1943 include 'COMMON.CONTROL'
1944 integer lenpre,lenpot,ilen,lentmp
1946 character*3 out1file_text,ucase
1949 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1950 call getenv_loc("PREFIX",prefix)
1952 call getenv_loc("POT",pot)
1953 call getenv_loc("DIRTMP",tmpdir)
1954 call getenv_loc("CURDIR",curdir)
1955 call getenv_loc("OUT1FILE",out1file_text)
1956 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1957 out1file_text=ucase(out1file_text)
1958 if (out1file_text(1:1).eq."Y") then
1961 out1file=fg_rank.gt.0
1966 if (lentmp.gt.0) then
1967 write (*,'(80(1h!))')
1968 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1969 write (*,'(80(1h!))')
1970 write (*,*)"All output files will be on node /tmp directory."
1972 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1973 if (me.eq.king) then
1974 write (*,*) "The master node is ",nodename
1975 else if (fg_rank.eq.0) then
1976 write (*,*) "I am the CG slave node ",nodename
1978 write (*,*) "I am the FG slave node ",nodename
1981 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1982 lenpre = lentmp+lenpre+1
1984 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1985 C Get the names and open the input files
1986 #if defined(WINIFL) || defined(WINPGI)
1987 open(1,file=pref_orig(:ilen(pref_orig))//
1988 & '.inp',status='old',readonly,shared)
1989 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1990 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1991 C Get parameter filenames and open the parameter files.
1992 call getenv_loc('BONDPAR',bondname)
1993 open (ibond,file=bondname,status='old',readonly,shared)
1994 call getenv_loc('THETPAR',thetname)
1995 open (ithep,file=thetname,status='old',readonly,shared)
1996 call getenv_loc('ROTPAR',rotname)
1997 open (irotam,file=rotname,status='old',readonly,shared)
1998 call getenv_loc('TORPAR',torname)
1999 open (itorp,file=torname,status='old',readonly,shared)
2000 call getenv_loc('TORDPAR',tordname)
2001 open (itordp,file=tordname,status='old',readonly,shared)
2002 call getenv_loc('FOURIER',fouriername)
2003 open (ifourier,file=fouriername,status='old',readonly,shared)
2004 call getenv_loc('ELEPAR',elename)
2005 open (ielep,file=elename,status='old',readonly,shared)
2006 call getenv_loc('SIDEPAR',sidename)
2007 open (isidep,file=sidename,status='old',readonly,shared)
2008 #elif (defined CRAY) || (defined AIX)
2009 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2011 c print *,"Processor",myrank," opened file 1"
2012 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2013 c print *,"Processor",myrank," opened file 9"
2014 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2015 C Get parameter filenames and open the parameter files.
2016 call getenv_loc('BONDPAR',bondname)
2017 open (ibond,file=bondname,status='old',action='read')
2018 c print *,"Processor",myrank," opened file IBOND"
2019 call getenv_loc('THETPAR',thetname)
2020 open (ithep,file=thetname,status='old',action='read')
2021 c print *,"Processor",myrank," opened file ITHEP"
2022 call getenv_loc('ROTPAR',rotname)
2023 open (irotam,file=rotname,status='old',action='read')
2024 c print *,"Processor",myrank," opened file IROTAM"
2025 call getenv_loc('TORPAR',torname)
2026 open (itorp,file=torname,status='old',action='read')
2027 c print *,"Processor",myrank," opened file ITORP"
2028 call getenv_loc('TORDPAR',tordname)
2029 open (itordp,file=tordname,status='old',action='read')
2030 c print *,"Processor",myrank," opened file ITORDP"
2031 call getenv_loc('SCCORPAR',sccorname)
2032 open (isccor,file=sccorname,status='old',action='read')
2033 c print *,"Processor",myrank," opened file ISCCOR"
2034 call getenv_loc('FOURIER',fouriername)
2035 open (ifourier,file=fouriername,status='old',action='read')
2036 c print *,"Processor",myrank," opened file IFOURIER"
2037 call getenv_loc('ELEPAR',elename)
2038 open (ielep,file=elename,status='old',action='read')
2039 c print *,"Processor",myrank," opened file IELEP"
2040 call getenv_loc('SIDEPAR',sidename)
2041 open (isidep,file=sidename,status='old',action='read')
2042 c print *,"Processor",myrank," opened file ISIDEP"
2043 c print *,"Processor",myrank," opened parameter files"
2045 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2046 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2047 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2048 C Get parameter filenames and open the parameter files.
2049 call getenv_loc('BONDPAR',bondname)
2050 open (ibond,file=bondname,status='old')
2051 call getenv_loc('THETPAR',thetname)
2052 open (ithep,file=thetname,status='old')
2053 call getenv_loc('ROTPAR',rotname)
2054 open (irotam,file=rotname,status='old')
2055 call getenv_loc('TORPAR',torname)
2056 open (itorp,file=torname,status='old')
2057 call getenv_loc('TORDPAR',tordname)
2058 open (itordp,file=tordname,status='old')
2059 call getenv_loc('SCCORPAR',sccorname)
2060 open (isccor,file=sccorname,status='old')
2061 call getenv_loc('FOURIER',fouriername)
2062 open (ifourier,file=fouriername,status='old')
2063 call getenv_loc('ELEPAR',elename)
2064 open (ielep,file=elename,status='old')
2065 call getenv_loc('SIDEPAR',sidename)
2066 open (isidep,file=sidename,status='old')
2068 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2070 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2071 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2072 C Get parameter filenames and open the parameter files.
2073 call getenv_loc('BONDPAR',bondname)
2074 open (ibond,file=bondname,status='old',readonly)
2075 call getenv_loc('THETPAR',thetname)
2076 open (ithep,file=thetname,status='old',readonly)
2077 call getenv_loc('ROTPAR',rotname)
2078 open (irotam,file=rotname,status='old',readonly)
2079 call getenv_loc('TORPAR',torname)
2080 open (itorp,file=torname,status='old',readonly)
2081 call getenv_loc('TORDPAR',tordname)
2082 open (itordp,file=tordname,status='old',readonly)
2083 call getenv_loc('SCCORPAR',sccorname)
2084 open (isccor,file=sccorname,status='old',readonly)
2086 call getenv_loc('THETPARPDB',thetname_pdb)
2087 print *,"thetname_pdb ",thetname_pdb
2088 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2089 print *,ithep_pdb," opened"
2091 call getenv_loc('FOURIER',fouriername)
2092 open (ifourier,file=fouriername,status='old',readonly)
2093 call getenv_loc('ELEPAR',elename)
2094 open (ielep,file=elename,status='old',readonly)
2095 call getenv_loc('SIDEPAR',sidename)
2096 open (isidep,file=sidename,status='old',readonly)
2097 call getenv_loc('LIPTRANPAR',liptranname)
2098 open (iliptranpar,file=liptranname,status='old',action='read')
2100 call getenv_loc('ROTPARPDB',rotname_pdb)
2101 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2106 C 8/9/01 In the newest version SCp interaction constants are read from a file
2107 C Use -DOLDSCP to use hard-coded constants instead.
2109 call getenv_loc('SCPPAR',scpname)
2110 #if defined(WINIFL) || defined(WINPGI)
2111 open (iscpp,file=scpname,status='old',readonly,shared)
2112 #elif (defined CRAY) || (defined AIX)
2113 open (iscpp,file=scpname,status='old',action='read')
2115 open (iscpp,file=scpname,status='old')
2117 open (iscpp,file=scpname,status='old',readonly)
2120 call getenv_loc('PATTERN',patname)
2121 #if defined(WINIFL) || defined(WINPGI)
2122 open (icbase,file=patname,status='old',readonly,shared)
2123 #elif (defined CRAY) || (defined AIX)
2124 open (icbase,file=patname,status='old',action='read')
2126 open (icbase,file=patname,status='old')
2128 open (icbase,file=patname,status='old',readonly)
2131 C Open output file only for CG processes
2132 c print *,"Processor",myrank," fg_rank",fg_rank
2133 if (fg_rank.eq.0) then
2135 if (nodes.eq.1) then
2138 npos = dlog10(dfloat(nodes-1))+1
2140 if (npos.lt.3) npos=3
2141 write (liczba,'(i1)') npos
2142 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2144 write (liczba,form) me
2145 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2146 & liczba(:ilen(liczba))
2147 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2149 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2151 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2152 & liczba(:ilen(liczba))//'.mol2'
2153 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2154 & liczba(:ilen(liczba))//'.stat'
2156 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2157 & //liczba(:ilen(liczba))//'.stat')
2158 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2161 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2162 & liczba(:ilen(liczba))//'.const'
2167 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2168 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2169 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2170 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2171 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2173 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2175 rest2name=prefix(:ilen(prefix))//'.rst'
2177 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2180 #if defined(AIX) || defined(PGI)
2181 if (me.eq.king .or. .not. out1file)
2182 & open(iout,file=outname,status='unknown')
2184 if (fg_rank.gt.0) then
2185 write (liczba,'(i3.3)') myrank/nfgtasks
2186 write (ll,'(bz,i3.3)') fg_rank
2187 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2192 open(igeom,file=intname,status='unknown',position='append')
2193 open(ipdb,file=pdbname,status='unknown')
2194 open(imol2,file=mol2name,status='unknown')
2195 open(istat,file=statname,status='unknown',position='append')
2197 c1out open(iout,file=outname,status='unknown')
2200 if (me.eq.king .or. .not.out1file)
2201 & open(iout,file=outname,status='unknown')
2203 if (fg_rank.gt.0) then
2204 write (liczba,'(i3.3)') myrank/nfgtasks
2205 write (ll,'(bz,i3.3)') fg_rank
2206 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2211 open(igeom,file=intname,status='unknown',access='append')
2212 open(ipdb,file=pdbname,status='unknown')
2213 open(imol2,file=mol2name,status='unknown')
2214 open(istat,file=statname,status='unknown',access='append')
2216 c1out open(iout,file=outname,status='unknown')
2219 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2220 csa_seed=prefix(:lenpre)//'.CSA.seed'
2221 csa_history=prefix(:lenpre)//'.CSA.history'
2222 csa_bank=prefix(:lenpre)//'.CSA.bank'
2223 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2224 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2225 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2226 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2227 csa_int=prefix(:lenpre)//'.int'
2228 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2229 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2230 csa_in=prefix(:lenpre)//'.CSA.in'
2231 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2234 write (iout,'(80(1h-))')
2235 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2236 write (iout,'(80(1h-))')
2237 write (iout,*) "Input file : ",
2238 & pref_orig(:ilen(pref_orig))//'.inp'
2239 write (iout,*) "Output file : ",
2240 & outname(:ilen(outname))
2242 write (iout,*) "Sidechain potential file : ",
2243 & sidename(:ilen(sidename))
2245 write (iout,*) "SCp potential file : ",
2246 & scpname(:ilen(scpname))
2248 write (iout,*) "Electrostatic potential file : ",
2249 & elename(:ilen(elename))
2250 write (iout,*) "Cumulant coefficient file : ",
2251 & fouriername(:ilen(fouriername))
2252 write (iout,*) "Torsional parameter file : ",
2253 & torname(:ilen(torname))
2254 write (iout,*) "Double torsional parameter file : ",
2255 & tordname(:ilen(tordname))
2256 write (iout,*) "SCCOR parameter file : ",
2257 & sccorname(:ilen(sccorname))
2258 write (iout,*) "Bond & inertia constant file : ",
2259 & bondname(:ilen(bondname))
2260 write (iout,*) "Bending parameter file : ",
2261 & thetname(:ilen(thetname))
2262 write (iout,*) "Rotamer parameter file : ",
2263 & rotname(:ilen(rotname))
2264 write (iout,*) "Thetpdb parameter file : ",
2265 & thetname_pdb(:ilen(thetname_pdb))
2266 write (iout,*) "Threading database : ",
2267 & patname(:ilen(patname))
2269 &write (iout,*)" DIRTMP : ",
2271 write (iout,'(80(1h-))')
2275 c----------------------------------------------------------------------------
2276 subroutine card_concat(card)
2277 implicit real*8 (a-h,o-z)
2278 include 'DIMENSIONS'
2279 include 'COMMON.IOUNITS'
2281 character*80 karta,ucase
2283 read (inp,'(a)') karta
2286 do while (karta(80:80).eq.'&')
2287 card=card(:ilen(card)+1)//karta(:79)
2288 read (inp,'(a)') karta
2291 card=card(:ilen(card)+1)//karta
2294 c----------------------------------------------------------------------------------
2296 implicit real*8 (a-h,o-z)
2297 include 'DIMENSIONS'
2298 include 'COMMON.CHAIN'
2299 include 'COMMON.IOUNITS'
2301 open(irest2,file=rest2name,status='unknown')
2302 read(irest2,*) totT,EK,potE,totE,t_bath
2305 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2308 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2311 read (irest2,*) iset
2316 c---------------------------------------------------------------------------------
2317 subroutine read_fragments
2318 implicit real*8 (a-h,o-z)
2319 include 'DIMENSIONS'
2323 include 'COMMON.SETUP'
2324 include 'COMMON.CHAIN'
2325 include 'COMMON.IOUNITS'
2327 include 'COMMON.CONTROL'
2328 read(inp,*) nset,nfrag,npair,nfrag_back
2329 if(me.eq.king.or..not.out1file)
2330 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2331 & " nfrag_back",nfrag_back
2333 read(inp,*) mset(iset)
2335 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2337 if(me.eq.king.or..not.out1file)
2338 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2339 & ifrag(2,i,iset), qinfrag(i,iset)
2342 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2344 if(me.eq.king.or..not.out1file)
2345 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2346 & ipair(2,i,iset), qinpair(i,iset)
2349 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2350 & wfrag_back(3,i,iset),
2351 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2352 if(me.eq.king.or..not.out1file)
2353 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2354 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2359 C---------------------------------------------------------------------------
2360 subroutine read_afminp
2361 implicit real*8 (a-h,o-z)
2362 include 'DIMENSIONS'
2366 include 'COMMON.SETUP'
2367 include 'COMMON.CONTROL'
2368 include 'COMMON.CHAIN'
2369 include 'COMMON.IOUNITS'
2370 include 'COMMON.SBRIDGE'
2371 character*320 afmcard
2373 call card_concat(afmcard)
2374 call readi(afmcard,"BEG",afmbeg,0)
2375 call readi(afmcard,"END",afmend,0)
2376 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2377 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2378 print *,'FORCE=' ,forceAFMconst
2379 CCCC NOW PROPERTIES FOR AFM
2382 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2384 distafminit=dsqrt(distafminit)
2385 print *,'initdist',distafminit
2389 c-------------------------------------------------------------------------------
2390 subroutine read_dist_constr
2391 implicit real*8 (a-h,o-z)
2392 include 'DIMENSIONS'
2396 include 'COMMON.SETUP'
2397 include 'COMMON.CONTROL'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2400 include 'COMMON.SBRIDGE'
2401 integer ifrag_(2,100),ipair_(2,100)
2402 double precision wfrag_(100),wpair_(100)
2403 character*500 controlcard
2405 write (iout,*) "Calling read_dist_constr"
2406 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2408 call card_concat(controlcard)
2409 call readi(controlcard,"NFRAG",nfrag_,0)
2410 call readi(controlcard,"NPAIR",npair_,0)
2411 call readi(controlcard,"NDIST",ndist_,0)
2412 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2413 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2414 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2415 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2416 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2417 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2418 c write (iout,*) "IFRAG"
2420 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2422 c write (iout,*) "IPAIR"
2424 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2428 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2429 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2430 & ifrag_(2,i)=nstart_sup+nsup-1
2431 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2433 if (wfrag_(i).gt.0.0d0) then
2434 do j=ifrag_(1,i),ifrag_(2,i)-1
2435 do k=j+1,ifrag_(2,i)
2436 c write (iout,*) "j",j," k",k
2438 if (constr_dist.eq.1) then
2443 forcon(nhpb)=wfrag_(i)
2444 else if (constr_dist.eq.2) then
2445 if (ddjk.le.dist_cut) then
2450 forcon(nhpb)=wfrag_(i)
2457 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2460 if (.not.out1file .or. me.eq.king)
2461 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2462 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2464 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2465 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2472 if (wpair_(i).gt.0.0d0) then
2480 do j=ifrag_(1,ii),ifrag_(2,ii)
2481 do k=ifrag_(1,jj),ifrag_(2,jj)
2485 forcon(nhpb)=wpair_(i)
2486 dhpb(nhpb)=dist(j,k)
2488 if (.not.out1file .or. me.eq.king)
2489 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2490 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2492 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2493 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2501 if (constr_dist.eq.11) then
2502 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2503 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2504 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2507 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2508 & ibecarb(i),forcon(nhpb+1)
2510 if (forcon(nhpb+1).gt.0.0d0) then
2512 if (ibecarb(i).gt.0) then
2513 ihpb(i)=ihpb(i)+nres
2514 jhpb(i)=jhpb(i)+nres
2516 if (dhpb(nhpb).eq.0.0d0)
2517 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2519 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2520 C if (forcon(nhpb+1).gt.0.0d0) then
2522 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2524 if (.not.out1file .or. me.eq.king)
2525 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2526 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2528 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2529 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2536 c-------------------------------------------------------------------------------
2538 subroutine flush(iu)
2543 subroutine flush(iu)
2548 c------------------------------------------------------------------------------
2549 subroutine copy_to_tmp(source)
2550 include "DIMENSIONS"
2551 include "COMMON.IOUNITS"
2552 character*(*) source
2553 character* 256 tmpfile
2557 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2558 inquire(file=tmpfile,exist=ex)
2560 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2561 & " to temporary directory..."
2562 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2563 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2567 c------------------------------------------------------------------------------
2568 subroutine move_from_tmp(source)
2569 include "DIMENSIONS"
2570 include "COMMON.IOUNITS"
2571 character*(*) source
2574 write (*,*) "Moving ",source(:ilen(source)),
2575 & " from temporary directory to working directory"
2576 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2577 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2580 c------------------------------------------------------------------------------
2581 subroutine random_init(seed)
2583 C Initialize random number generator
2585 implicit real*8 (a-h,o-z)
2586 include 'DIMENSIONS'
2589 logical OKRandom, prng_restart
2591 integer iseed_array(4)
2593 include 'COMMON.IOUNITS'
2594 include 'COMMON.TIME1'
2595 include 'COMMON.THREAD'
2596 include 'COMMON.SBRIDGE'
2597 include 'COMMON.CONTROL'
2598 include 'COMMON.MCM'
2599 include 'COMMON.MAP'
2600 include 'COMMON.HEADER'
2601 include 'COMMON.CSA'
2602 include 'COMMON.CHAIN'
2603 include 'COMMON.MUCA'
2605 include 'COMMON.FFIELD'
2606 include 'COMMON.SETUP'
2607 iseed=-dint(dabs(seed))
2608 if (iseed.eq.0) then
2609 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2610 & 'Random seed undefined. The program will stop.'
2611 write (*,'(/80(1h*)/20x,a/80(1h*))')
2612 & 'Random seed undefined. The program will stop.'
2614 call mpi_finalize(mpi_comm_world,ierr)
2616 stop 'Bad random seed.'
2619 if (fg_rank.eq.0) then
2623 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2624 OKRandom = prng_restart(me,iseed)
2627 tmp=65536.0d0**(4-i)
2628 iseed_array(i) = dint(seed/tmp)
2629 seed=seed-iseed_array(i)*tmp
2632 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2633 & (iseed_array(i),i=1,4)
2634 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2635 & (iseed_array(i),i=1,4)
2636 OKRandom = prng_restart(me,iseed_array)
2639 c r1 = prng_next(me)
2640 r1=ran_number(0.0D0,1.0D0)
2642 & write (iout,*) 'ran_num',r1
2643 if (r1.lt.0.0d0) OKRandom=.false.
2645 if (.not.OKRandom) then
2646 write (iout,*) 'PRNG IS NOT WORKING!!!'
2647 print *,'PRNG IS NOT WORKING!!!'
2650 call mpi_abort(mpi_comm_world,error_msg,ierr)
2653 write (iout,*) 'too many processors for parallel prng'
2654 write (*,*) 'too many processors for parallel prng'
2662 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)