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,2800.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,'TUBEMOD',tubelog,0)
149 write (iout,*) TUBElog,"TUBEMODE"
150 if (TUBElog.gt.0) then
151 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
152 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
153 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
155 call readi(controlcard,'IPRINT',iprint,0)
156 C SHIELD keyword sets if the shielding effect of side-chains is used
157 C 0 denots no shielding is used all peptide are equally despite the
158 C solvent accesible area
159 C 1 the newly introduced function
160 C 2 reseved for further possible developement
161 call readi(controlcard,'SHIELD',shield_mode,0)
162 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
163 write(iout,*) "shield_mode",shield_mode
165 call readi(controlcard,'TORMODE',tor_mode,0)
166 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
167 write(iout,*) "torsional and valence angle mode",tor_mode
168 call readi(controlcard,'MAXGEN',maxgen,10000)
169 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
170 call readi(controlcard,"KDIAG",kdiag,0)
171 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
172 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
173 & write (iout,*) "RESCALE_MODE",rescale_mode
174 split_ene=index(controlcard,'SPLIT_ENE').gt.0
175 if (index(controlcard,'REGULAR').gt.0.0D0) then
176 call reada(controlcard,'WEIDIS',weidis,0.1D0)
180 if (index(controlcard,'CHECKGRAD').gt.0) then
182 if (index(controlcard,'CART').gt.0) then
184 elseif (index(controlcard,'CARINT').gt.0) then
189 call reada(controlcard,'DELTA',aincr,1.0d-7)
190 elseif (index(controlcard,'THREAD').gt.0) then
192 call readi(controlcard,'THREAD',nthread,0)
193 if (nthread.gt.0) then
194 call reada(controlcard,'WEIDIS',weidis,0.1D0)
197 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
198 stop 'Error termination in Read_Control.'
200 else if (index(controlcard,'MCMA').gt.0) then
202 else if (index(controlcard,'MCEE').gt.0) then
204 else if (index(controlcard,'MULTCONF').gt.0) then
206 else if (index(controlcard,'MAP').gt.0) then
208 call readi(controlcard,'MAP',nmap,0)
209 else if (index(controlcard,'CSA').gt.0) then
211 crc else if (index(controlcard,'ZSCORE').gt.0) then
213 crc ZSCORE is rm from UNRES, modecalc=9 is available
216 cfcm else if (index(controlcard,'MCMF').gt.0) then
218 else if (index(controlcard,'SOFTREG').gt.0) then
220 else if (index(controlcard,'CHECK_BOND').gt.0) then
222 else if (index(controlcard,'TEST').gt.0) then
224 else if (index(controlcard,'MD').gt.0) then
226 else if (index(controlcard,'RE ').gt.0) then
230 lmuca=index(controlcard,'MUCA').gt.0
231 call readi(controlcard,'MUCADYN',mucadyn,0)
232 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
233 if (lmuca .and. (me.eq.king .or. .not.out1file ))
235 write (iout,*) 'MUCADYN=',mucadyn
236 write (iout,*) 'MUCASMOOTH=',muca_smooth
239 iscode=index(controlcard,'ONE_LETTER')
240 indphi=index(controlcard,'PHI')
241 indback=index(controlcard,'BACK')
242 iranconf=index(controlcard,'RAND_CONF')
243 i2ndstr=index(controlcard,'USE_SEC_PRED')
244 gradout=index(controlcard,'GRADOUT').gt.0
245 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
246 C DISTCHAINMAX become obsolete for periodic boundry condition
247 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
248 C Reading the dimensions of box in x,y,z coordinates
249 call reada(controlcard,'BOXX',boxxsize,100.0d0)
250 call reada(controlcard,'BOXY',boxysize,100.0d0)
251 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
252 c Cutoff range for interactions
253 call reada(controlcard,"R_CUT",r_cut,15.0d0)
254 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
255 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
256 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
257 if (lipthick.gt.0.0d0) then
258 bordliptop=(boxzsize+lipthick)/2.0
259 bordlipbot=bordliptop-lipthick
261 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
262 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
263 buflipbot=bordlipbot+lipbufthick
264 bufliptop=bordliptop-lipbufthick
265 if ((lipbufthick*2.0d0).gt.lipthick)
266 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
268 write(iout,*) "bordliptop=",bordliptop
269 write(iout,*) "bordlipbot=",bordlipbot
270 write(iout,*) "bufliptop=",bufliptop
271 write(iout,*) "buflipbot=",buflipbot
272 write (iout,*) "SHIELD MODE",shield_mode
273 if (shield_mode.gt.0) then
275 C VSolvSphere the volume of solving sphere
277 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
278 C there will be no distinction between proline peptide group and normal peptide
279 C group in case of shielding parameters
280 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
281 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
282 write (iout,*) VSolvSphere,VSolvSphere_div
283 C long axis of side chain
285 long_r_sidechain(i)=vbldsc0(1,i)
286 short_r_sidechain(i)=sigma0(i)
290 if (me.eq.king .or. .not.out1file )
291 & write (iout,*) "DISTCHAINMAX",distchainmax
293 if(me.eq.king.or..not.out1file)
294 & write (iout,'(2a)') diagmeth(kdiag),
295 & ' routine used to diagonalize matrices.'
298 c--------------------------------------------------------------------------
299 subroutine read_REMDpar
303 implicit real*8 (a-h,o-z)
305 include 'COMMON.IOUNITS'
306 include 'COMMON.TIME1'
309 include 'COMMON.LANGEVIN'
311 include 'COMMON.LANGEVIN.lang0'
313 include 'COMMON.INTERACT'
314 include 'COMMON.NAMES'
316 include 'COMMON.REMD'
317 include 'COMMON.CONTROL'
318 include 'COMMON.SETUP'
320 character*320 controlcard
321 character*3200 controlcard1
322 integer iremd_m_total
324 if(me.eq.king.or..not.out1file)
325 & write (iout,*) "REMD setup"
327 call card_concat(controlcard)
328 call readi(controlcard,"NREP",nrep,3)
329 call readi(controlcard,"NSTEX",nstex,1000)
330 call reada(controlcard,"RETMIN",retmin,10.0d0)
331 call reada(controlcard,"RETMAX",retmax,1000.0d0)
332 mremdsync=(index(controlcard,'SYNC').gt.0)
333 call readi(controlcard,"NSYN",i_sync_step,100)
334 restart1file=(index(controlcard,'REST1FILE').gt.0)
335 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
336 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
337 if(max_cache_traj_use.gt.max_cache_traj)
338 & max_cache_traj_use=max_cache_traj
339 if(me.eq.king.or..not.out1file) then
340 cd if (traj1file) then
341 crc caching is in testing - NTWX is not ignored
342 cd write (iout,*) "NTWX value is ignored"
343 cd write (iout,*) " trajectory is stored to one file by master"
344 cd write (iout,*) " before exchange at NSTEX intervals"
346 write (iout,*) "NREP= ",nrep
347 write (iout,*) "NSTEX= ",nstex
348 write (iout,*) "SYNC= ",mremdsync
349 write (iout,*) "NSYN= ",i_sync_step
350 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
353 if (index(controlcard,'TLIST').gt.0) then
355 call card_concat(controlcard1)
356 read(controlcard1,*) (remd_t(i),i=1,nrep)
357 if(me.eq.king.or..not.out1file)
358 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
361 if (index(controlcard,'MLIST').gt.0) then
363 call card_concat(controlcard1)
364 read(controlcard1,*) (remd_m(i),i=1,nrep)
365 if(me.eq.king.or..not.out1file) then
366 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
369 iremd_m_total=iremd_m_total+remd_m(i)
371 write (iout,*) 'Total number of replicas ',iremd_m_total
374 if(me.eq.king.or..not.out1file)
375 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
378 c--------------------------------------------------------------------------
379 subroutine read_MDpar
383 implicit real*8 (a-h,o-z)
385 include 'COMMON.IOUNITS'
386 include 'COMMON.TIME1'
389 include 'COMMON.LANGEVIN'
391 include 'COMMON.LANGEVIN.lang0'
393 include 'COMMON.INTERACT'
394 include 'COMMON.NAMES'
396 include 'COMMON.SETUP'
397 include 'COMMON.CONTROL'
398 include 'COMMON.SPLITELE'
400 character*320 controlcard
402 call card_concat(controlcard)
403 call readi(controlcard,"NSTEP",n_timestep,1000000)
404 call readi(controlcard,"NTWE",ntwe,100)
405 call readi(controlcard,"NTWX",ntwx,1000)
406 call reada(controlcard,"DT",d_time,1.0d-1)
407 call reada(controlcard,"DVMAX",dvmax,2.0d1)
408 call reada(controlcard,"DAMAX",damax,1.0d1)
409 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
410 call readi(controlcard,"LANG",lang,0)
411 RESPA = index(controlcard,"RESPA") .gt. 0
412 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
413 ntime_split0=ntime_split
414 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
415 ntime_split0=ntime_split
416 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
417 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
418 rest = index(controlcard,"REST").gt.0
419 tbf = index(controlcard,"TBF").gt.0
420 usampl = index(controlcard,"USAMPL").gt.0
422 mdpdb = index(controlcard,"MDPDB").gt.0
423 call reada(controlcard,"T_BATH",t_bath,300.0d0)
424 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
425 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
426 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
427 if (count_reset_moment.eq.0) count_reset_moment=1000000000
428 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
429 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
430 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
431 if (count_reset_vel.eq.0) count_reset_vel=1000000000
432 large = index(controlcard,"LARGE").gt.0
433 print_compon = index(controlcard,"PRINT_COMPON").gt.0
434 rattle = index(controlcard,"RATTLE").gt.0
435 c if performing umbrella sampling, fragments constrained are read from the fragment file
441 if(me.eq.king.or..not.out1file) then
443 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
445 write (iout,'(a)') "The units are:"
446 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
447 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
448 & " acceleration: angstrom/(48.9 fs)**2"
449 write (iout,'(a)') "energy: kcal/mol, temperature: K"
451 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
452 write (iout,'(a60,f10.5,a)')
453 & "Initial time step of numerical integration:",d_time,
455 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
457 write (iout,'(2a,i4,a)')
458 & "A-MTS algorithm used; initial time step for fast-varying",
459 & " short-range forces split into",ntime_split," steps."
460 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
461 & r_cut," lambda",rlamb
463 write (iout,'(2a,f10.5)')
464 & "Maximum acceleration threshold to reduce the time step",
465 & "/increase split number:",damax
466 write (iout,'(2a,f10.5)')
467 & "Maximum predicted energy drift to reduce the timestep",
468 & "/increase split number:",edriftmax
469 write (iout,'(a60,f10.5)')
470 & "Maximum velocity threshold to reduce velocities:",dvmax
471 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
472 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
473 if (rattle) write (iout,'(a60)')
474 & "Rattle algorithm used to constrain the virtual bonds"
478 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
479 call reada(controlcard,"RWAT",rwat,1.4d0)
480 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
481 surfarea=index(controlcard,"SURFAREA").gt.0
482 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
483 if(me.eq.king.or..not.out1file)then
484 write (iout,'(/a,$)') "Langevin dynamics calculation"
487 & " with direct integration of Langevin equations"
488 else if (lang.eq.2) then
489 write (iout,'(a/)') " with TINKER stochasic MD integrator"
490 else if (lang.eq.3) then
491 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
492 else if (lang.eq.4) then
493 write (iout,'(a/)') " in overdamped mode"
495 write (iout,'(//a,i5)')
496 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
499 write (iout,'(a60,f10.5)') "Temperature:",t_bath
500 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
501 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
502 write (iout,'(a60,f10.5)')
503 & "Scaling factor of the friction forces:",scal_fric
504 if (surfarea) write (iout,'(2a,i10,a)')
505 & "Friction coefficients will be scaled by solvent-accessible",
506 & " surface area every",reset_fricmat," steps."
508 c Calculate friction coefficients and bounds of stochastic forces
509 eta=6*pi*cPoise*etawat
510 if(me.eq.king.or..not.out1file)
511 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
513 gamp=scal_fric*(pstok+rwat)*eta
514 stdfp=dsqrt(2*Rb*t_bath/d_time)
516 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
517 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
519 if(me.eq.king.or..not.out1file)then
520 write (iout,'(/2a/)')
521 & "Radii of site types and friction coefficients and std's of",
522 & " stochastic forces of fully exposed sites"
523 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
525 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
526 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
530 if(me.eq.king.or..not.out1file)then
531 write (iout,'(a)') "Berendsen bath calculation"
532 write (iout,'(a60,f10.5)') "Temperature:",t_bath
533 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
535 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
536 & count_reset_moment," steps"
538 & write (iout,'(a,i10,a)')
539 & "Velocities will be reset at random every",count_reset_vel,
543 if(me.eq.king.or..not.out1file)
544 & write (iout,'(a31)') "Microcanonical mode calculation"
546 if(me.eq.king.or..not.out1file)then
547 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
549 write(iout,*) "MD running with constraints."
550 write(iout,*) "Equilibration time ", eq_time, " mtus."
551 write(iout,*) "Constraining ", nfrag," fragments."
552 write(iout,*) "Length of each fragment, weight and q0:"
554 write (iout,*) "Set of restraints #",iset
556 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
557 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
559 write(iout,*) "constraints between ", npair, "fragments."
560 write(iout,*) "constraint pairs, weights and q0:"
562 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
563 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
565 write(iout,*) "angle constraints within ", nfrag_back,
566 & "backbone fragments."
567 write(iout,*) "fragment, weights:"
569 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
570 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
571 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
574 iset=mod(kolor,nset)+1
577 if(me.eq.king.or..not.out1file)
578 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
581 c------------------------------------------------------------------------------
584 C Read molecular data.
586 implicit real*8 (a-h,o-z)
592 include 'COMMON.IOUNITS'
595 include 'COMMON.INTERACT'
596 include 'COMMON.LOCAL'
597 include 'COMMON.NAMES'
598 include 'COMMON.CHAIN'
599 include 'COMMON.FFIELD'
600 include 'COMMON.SBRIDGE'
601 include 'COMMON.HEADER'
602 include 'COMMON.CONTROL'
603 include 'COMMON.DBASE'
604 include 'COMMON.THREAD'
605 include 'COMMON.CONTACTS'
606 include 'COMMON.TORCNSTR'
607 include 'COMMON.TIME1'
608 include 'COMMON.BOUNDS'
610 include 'COMMON.SETUP'
611 include 'COMMON.SHIELD'
612 character*4 sequence(maxres)
614 double precision x(maxvar)
615 character*256 pdbfile
616 character*400 weightcard
617 character*80 weightcard_t,ucase
618 dimension itype_pdb(maxres)
619 common /pizda/ itype_pdb
620 logical seq_comp,fail
621 double precision energia(0:n_ene)
627 C Read weights of the subsequent energy terms.
628 call card_concat(weightcard)
629 call reada(weightcard,'WLONG',wlong,1.0D0)
630 call reada(weightcard,'WSC',wsc,wlong)
631 call reada(weightcard,'WSCP',wscp,wlong)
632 call reada(weightcard,'WELEC',welec,1.0D0)
633 call reada(weightcard,'WVDWPP',wvdwpp,welec)
634 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
635 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
636 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
637 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
638 call reada(weightcard,'WTURN3',wturn3,1.0D0)
639 call reada(weightcard,'WTURN4',wturn4,1.0D0)
640 call reada(weightcard,'WTURN6',wturn6,1.0D0)
641 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
642 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
643 call reada(weightcard,'WBOND',wbond,1.0D0)
644 call reada(weightcard,'WTOR',wtor,1.0D0)
645 call reada(weightcard,'WTORD',wtor_d,1.0D0)
646 call reada(weightcard,'WANG',wang,1.0D0)
647 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
648 call reada(weightcard,'SCAL14',scal14,0.4D0)
649 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
650 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
651 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
652 call reada(weightcard,'TEMP0',temp0,300.0d0)
653 call reada(weightcard,'WSHIELD',wshield,1.0d0)
654 call reada(weightcard,'WLT',wliptran,0.0D0)
655 call reada(weightcard,'WTUBE',wtube,1.0D0)
656 if (index(weightcard,'SOFT').gt.0) ipot=6
657 C 12/1/95 Added weight for the multi-body term WCORR
658 call reada(weightcard,'WCORRH',wcorr,1.0D0)
659 if (wcorr4.gt.0.0d0) wcorr=wcorr4
679 if(me.eq.king.or..not.out1file)
680 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
681 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
683 10 format (/'Energy-term weights (unscaled):'//
684 & 'WSCC= ',f10.6,' (SC-SC)'/
685 & 'WSCP= ',f10.6,' (SC-p)'/
686 & 'WELEC= ',f10.6,' (p-p electr)'/
687 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
688 & 'WBOND= ',f10.6,' (stretching)'/
689 & 'WANG= ',f10.6,' (bending)'/
690 & 'WSCLOC= ',f10.6,' (SC local)'/
691 & 'WTOR= ',f10.6,' (torsional)'/
692 & 'WTORD= ',f10.6,' (double torsional)'/
693 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
694 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
695 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
696 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
697 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
698 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
699 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
700 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
701 & 'WTURN6= ',f10.6,' (turns, 6th order)')
702 if(me.eq.king.or..not.out1file)then
703 if (wcorr4.gt.0.0d0) then
704 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
705 & 'between contact pairs of peptide groups'
706 write (iout,'(2(a,f5.3/))')
707 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
708 & 'Range of quenching the correlation terms:',2*delt_corr
709 else if (wcorr.gt.0.0d0) then
710 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
711 & 'between contact pairs of peptide groups'
713 write (iout,'(a,f8.3)')
714 & 'Scaling factor of 1,4 SC-p interactions:',scal14
715 write (iout,'(a,f8.3)')
716 & 'General scaling factor of SC-p interactions:',scalscp
718 r0_corr=cutoff_corr-delt_corr
720 aad(i,1)=scalscp*aad(i,1)
721 aad(i,2)=scalscp*aad(i,2)
722 bad(i,1)=scalscp*bad(i,1)
723 bad(i,2)=scalscp*bad(i,2)
725 call rescale_weights(t_bath)
726 if(me.eq.king.or..not.out1file)
727 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
728 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
730 22 format (/'Energy-term weights (scaled):'//
731 & 'WSCC= ',f10.6,' (SC-SC)'/
732 & 'WSCP= ',f10.6,' (SC-p)'/
733 & 'WELEC= ',f10.6,' (p-p electr)'/
734 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
735 & 'WBOND= ',f10.6,' (stretching)'/
736 & 'WANG= ',f10.6,' (bending)'/
737 & 'WSCLOC= ',f10.6,' (SC local)'/
738 & 'WTOR= ',f10.6,' (torsional)'/
739 & 'WTORD= ',f10.6,' (double torsional)'/
740 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
741 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
742 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
743 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
744 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
745 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
746 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
747 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
748 & 'WTURN6= ',f10.6,' (turns, 6th order)')
749 if(me.eq.king.or..not.out1file)
750 & write (iout,*) "Reference temperature for weights calculation:",
752 call reada(weightcard,"D0CM",d0cm,3.78d0)
753 call reada(weightcard,"AKCM",akcm,15.1d0)
754 call reada(weightcard,"AKTH",akth,11.0d0)
755 call reada(weightcard,"AKCT",akct,12.0d0)
756 call reada(weightcard,"V1SS",v1ss,-1.08d0)
757 call reada(weightcard,"V2SS",v2ss,7.61d0)
758 call reada(weightcard,"V3SS",v3ss,13.7d0)
759 call reada(weightcard,"EBR",ebr,-5.50D0)
760 call reada(weightcard,"ATRISS",atriss,0.301D0)
761 call reada(weightcard,"BTRISS",btriss,0.021D0)
762 call reada(weightcard,"CTRISS",ctriss,1.001D0)
763 call reada(weightcard,"DTRISS",dtriss,1.001D0)
764 write (iout,*) "ATRISS=", atriss
765 write (iout,*) "BTRISS=", btriss
766 write (iout,*) "CTRISS=", ctriss
767 write (iout,*) "DTRISS=", dtriss
768 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
770 dyn_ss_mask(i)=.false.
774 dyn_ssbond_ij(i,j)=1.0d300
777 call reada(weightcard,"HT",Ht,0.0D0)
779 ss_depth=ebr/wsc-0.25*eps(1,1)
780 Ht=Ht/wsc-0.25*eps(1,1)
781 akcm=akcm*wstrain/wsc
782 akth=akth*wstrain/wsc
783 akct=akct*wstrain/wsc
784 v1ss=v1ss*wstrain/wsc
785 v2ss=v2ss*wstrain/wsc
786 v3ss=v3ss*wstrain/wsc
788 if (wstrain.ne.0.0) then
789 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
794 write (iout,*) "wshield,", wshield
795 if(me.eq.king.or..not.out1file) then
796 write (iout,*) "Parameters of the SS-bond potential:"
797 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
799 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
800 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
801 write (iout,*)" HT",Ht
802 print *,'indpdb=',indpdb,' pdbref=',pdbref
804 if (indpdb.gt.0 .or. pdbref) then
805 read(inp,'(a)') pdbfile
806 if(me.eq.king.or..not.out1file)
807 & write (iout,'(2a)') 'PDB data will be read from file ',
808 & pdbfile(:ilen(pdbfile))
809 open(ipdbin,file=pdbfile,status='old',err=33)
811 33 write (iout,'(a)') 'Error opening PDB file.'
814 c write (iout,*) 'Begin reading pdb data'
817 c write (iout,*) 'Finished reading pdb data'
819 if(me.eq.king.or..not.out1file)
820 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
821 & ' nstart_sup=',nstart_sup
823 itype_pdb(i)=itype(i)
827 nct=nstart_sup+nsup-1
828 call contact(.false.,ncont_ref,icont_ref,co)
831 if(me.eq.king.or..not.out1file)
832 & write(iout,*)'Adding sidechains'
836 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
839 do while (fail.and.nsi.le.maxsi)
840 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
843 if(fail) write(iout,*)'Adding sidechain failed for res ',
844 & i,' after ',nsi,' trials'
849 if (indpdb.eq.0) then
850 C Read sequence if not taken from the pdb file.
852 c print *,'nres=',nres
853 if (iscode.gt.0) then
854 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
856 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
858 C Convert sequence to numeric code
860 itype(i)=rescode(i,sequence(i),iscode)
862 C Assign initial virtual bond lengths
868 vbld(i+nres)=dsc(iabs(itype(i)))
869 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
870 c write (iout,*) "i",i," itype",itype(i),
871 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
875 c print '(20i4)',(itype(i),i=1,nres)
878 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
880 if (itype(i).eq.ntyp1) then
884 else if (iabs(itype(i+1)).ne.20) then
886 else if (iabs(itype(i)).ne.20) then
893 if(me.eq.king.or..not.out1file)then
894 write (iout,*) "ITEL"
896 write (iout,*) i,itype(i),itel(i)
898 print *,'Call Read_Bridge.'
901 C 8/13/98 Set limits to generating the dihedral angles
906 read (inp,*) ndih_constr
907 if (ndih_constr.gt.0) then
909 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
911 if(me.eq.king.or..not.out1file)then
913 & 'There are',ndih_constr,' constraints on phi angles.'
915 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
920 phi0(i)=deg2rad*phi0(i)
921 drange(i)=deg2rad*drange(i)
923 C if(me.eq.king.or..not.out1file)
924 C & write (iout,*) 'FTORS',ftors
927 phibound(1,ii) = phi0(i)-drange(i)
928 phibound(2,ii) = phi0(i)+drange(i)
931 C first setting the theta boundaries to 0 to pi
932 C this mean that there is no energy penalty for any angle occuring this can be applied
933 C for generate random conformation but is not implemented in this way
938 C begin reading theta constrains this is quartic constrains allowing to
939 C have smooth second derivative
940 if (with_theta_constr) then
941 C with_theta_constr is keyword allowing for occurance of theta constrains
942 read (inp,*) ntheta_constr
943 C ntheta_constr is the number of theta constrains
944 if (ntheta_constr.gt.0) then
946 read (inp,*) (itheta_constr(i),theta_constr0(i),
947 & theta_drange(i),for_thet_constr(i),
949 C the above code reads from 1 to ntheta_constr
950 C itheta_constr(i) residue i for which is theta_constr
951 C theta_constr0 the global minimum value
952 C theta_drange is range for which there is no energy penalty
953 C for_thet_constr is the force constant for quartic energy penalty
955 if(me.eq.king.or..not.out1file)then
957 & 'There are',ntheta_constr,' constraints on phi angles.'
959 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
965 theta_constr0(i)=deg2rad*theta_constr0(i)
966 theta_drange(i)=deg2rad*theta_drange(i)
968 C if(me.eq.king.or..not.out1file)
969 C & write (iout,*) 'FTORS',ftors
970 C do i=1,ntheta_constr
971 C ii = itheta_constr(i)
972 C thetabound(1,ii) = phi0(i)-drange(i)
973 C thetabound(2,ii) = phi0(i)+drange(i)
975 endif ! ntheta_constr.gt.0
976 endif! with_theta_constr
978 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
979 C write (iout,*) "with_dihed_constr ",with_dihed_constr
984 write (iout,'(a)') 'Boundaries in phi angle sampling:'
986 write (iout,'(a3,i5,2f10.1)')
987 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
993 cd print *,'NNT=',NNT,' NCT=',NCT
994 if (itype(1).eq.ntyp1) nnt=2
995 if (itype(nres).eq.ntyp1) nct=nct-1
997 if(me.eq.king.or..not.out1file)
998 & write (iout,'(a,i3)') 'nsup=',nsup
1000 if (nsup.le.(nct-nnt+1)) then
1001 do i=0,nct-nnt+1-nsup
1002 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1008 & 'Error - sequences to be superposed do not match.'
1011 do i=0,nsup-(nct-nnt+1)
1012 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1014 nstart_sup=nstart_sup+i
1020 & 'Error - sequences to be superposed do not match.'
1023 if (nsup.eq.0) nsup=nct-nnt
1024 if (nstart_sup.eq.0) nstart_sup=nnt
1025 if (nstart_seq.eq.0) nstart_seq=nnt
1026 if(me.eq.king.or..not.out1file)
1027 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1028 & ' nstart_seq=',nstart_seq
1030 c--- Zscore rms -------
1031 if (nz_start.eq.0) nz_start=nnt
1032 if (nz_end.eq.0 .and. nsup.gt.0) then
1034 else if (nz_end.eq.0) then
1037 if(me.eq.king.or..not.out1file)then
1038 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1039 write (iout,*) 'IZ_SC=',iz_sc
1041 c----------------------
1044 if (.not.pdbref) then
1045 call read_angles(inp,*38)
1047 38 write (iout,'(a)') 'Error reading reference structure.'
1049 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1050 stop 'Error reading reference structure'
1052 39 call chainbuild_extconf
1054 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1061 cref(j,i,kkk)=c(j,i)
1064 call contact(.true.,ncont_ref,icont_ref,co)
1068 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1070 if (constr_dist.gt.0) call read_dist_constr
1071 write (iout,*) "After read_dist_constr nhpb",nhpb
1072 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1074 if(me.eq.king.or..not.out1file)
1075 & write (iout,*) 'Contact order:',co
1077 if(me.eq.king.or..not.out1file)
1078 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1081 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1083 if(me.eq.king.or..not.out1file)
1084 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1085 & icont_ref(1,i),' ',
1086 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1090 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1091 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1092 & modecalc.ne.10) then
1093 C If input structure hasn't been supplied from the PDB file read or generate
1095 if (iranconf.eq.0 .and. .not. extconf) then
1096 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1097 & write (iout,'(a)') 'Initial geometry will be read in.'
1099 read(inp,'(8f10.5)',end=36,err=36)
1100 & ((c(l,k),l=1,3),k=1,nres),
1101 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1102 write (iout,*) "Exit READ_CART"
1103 write (iout,'(8f10.5)')
1104 & ((c(l,k),l=1,3),k=1,nres),
1105 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1106 call int_from_cart1(.true.)
1107 write (iout,*) "Finish INT_TO_CART"
1110 dc(j,i)=c(j,i+1)-c(j,i)
1111 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1115 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1117 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1118 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1124 call read_angles(inp,*36)
1125 call chainbuild_extconf
1128 36 write (iout,'(a)') 'Error reading angle file.'
1130 call mpi_finalize( MPI_COMM_WORLD,IERR )
1132 stop 'Error reading angle file.'
1134 else if (extconf) then
1135 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1136 & write (iout,'(a)') 'Extended chain initial geometry.'
1138 theta(i)=90d0*deg2rad
1141 phi(i)=180d0*deg2rad
1144 alph(i)=110d0*deg2rad
1147 omeg(i)=-120d0*deg2rad
1148 if (itype(i).le.0) omeg(i)=-omeg(i)
1150 call chainbuild_extconf
1152 if(me.eq.king.or..not.out1file)
1153 & write (iout,'(a)') 'Random-generated initial geometry.'
1157 if (me.eq.king .or. fg_rank.eq.0 .and. (
1158 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1162 call gen_rand_conf(itmp,*30)
1164 30 write (iout,*) 'Failed to generate random conformation',
1165 & ', itrial=',itrial
1166 write (*,*) 'Processor:',me,
1167 & ' Failed to generate random conformation',
1177 write (iout,'(a,i3,a)') 'Processor:',me,
1178 & ' error in generating random conformation.'
1179 write (*,'(a,i3,a)') 'Processor:',me,
1180 & ' error in generating random conformation.'
1183 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1188 & ' error in generating random conformation.'
1193 elseif (modecalc.eq.4) then
1194 read (inp,'(a)') intinname
1195 open (intin,file=intinname,status='old',err=333)
1196 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1197 & write (iout,'(a)') 'intinname',intinname
1198 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1200 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1202 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1204 stop 'Error opening angle file.'
1208 C Generate distance constraints, if the PDB structure is to be regularized.
1209 if (nthread.gt.0) then
1210 call read_threadbase
1213 if (me.eq.king .or. .not. out1file)
1215 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1216 write (iout,'(/a,i3,a)')
1217 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1218 write (iout,'(20i4)') (iss(i),i=1,ns)
1220 write(iout,*)"Running with dynamic disulfide-bond formation"
1222 write (iout,'(/a/)') 'Pre-formed links are:'
1228 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1229 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1235 if (ns.gt.0.and.dyn_ss) then
1239 forcon(i-nss)=forcon(i)
1246 dyn_ss_mask(iss(i))=.true.
1249 if (i2ndstr.gt.0) call secstrp2dihc
1250 c call geom_to_var(nvar,x)
1251 c call etotal(energia(0))
1252 c call enerprint(energia(0))
1253 c call briefout(0,etot)
1255 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1256 cd write (iout,'(a)') 'Variable list:'
1257 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1259 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1260 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1261 & 'Processor',myrank,': end reading molecular data.'
1266 c--------------------------------------------------------------------------
1267 logical function seq_comp(itypea,itypeb,length)
1269 integer length,itypea(length),itypeb(length)
1272 if (itypea(i).ne.itypeb(i)) then
1280 c-----------------------------------------------------------------------------
1281 subroutine read_bridge
1282 C Read information about disulfide bridges.
1283 implicit real*8 (a-h,o-z)
1284 include 'DIMENSIONS'
1288 include 'COMMON.IOUNITS'
1289 include 'COMMON.GEO'
1290 include 'COMMON.VAR'
1291 include 'COMMON.INTERACT'
1292 include 'COMMON.LOCAL'
1293 include 'COMMON.NAMES'
1294 include 'COMMON.CHAIN'
1295 include 'COMMON.FFIELD'
1296 include 'COMMON.SBRIDGE'
1297 include 'COMMON.HEADER'
1298 include 'COMMON.CONTROL'
1299 include 'COMMON.DBASE'
1300 include 'COMMON.THREAD'
1301 include 'COMMON.TIME1'
1302 include 'COMMON.SETUP'
1303 C Read bridging residues.
1304 read (inp,*) ns,(iss(i),i=1,ns)
1306 if(me.eq.king.or..not.out1file)
1307 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1308 C Check whether the specified bridging residues are cystines.
1310 if (itype(iss(i)).ne.1) then
1311 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1312 & 'Do you REALLY think that the residue ',
1313 & restyp(itype(iss(i))),i,
1314 & ' can form a disulfide bridge?!!!'
1315 write (*,'(2a,i3,a)')
1316 & 'Do you REALLY think that the residue ',
1317 & restyp(itype(iss(i))),i,
1318 & ' can form a disulfide bridge?!!!'
1320 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1325 C Read preformed bridges.
1327 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1329 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1332 C Check if the residues involved in bridges are in the specified list of
1333 C bridging residues.
1336 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1337 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1338 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1339 & ' contains residues present in other pairs.'
1340 write (*,'(a,i3,a)') 'Disulfide pair',i,
1341 & ' contains residues present in other pairs.'
1343 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1349 if (ihpb(i).eq.iss(j)) goto 10
1351 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1354 if (jhpb(i).eq.iss(j)) goto 20
1356 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1362 ihpb(i)=ihpb(i)+nres
1363 jhpb(i)=jhpb(i)+nres
1369 c----------------------------------------------------------------------------
1370 subroutine read_x(kanal,*)
1371 implicit real*8 (a-h,o-z)
1372 include 'DIMENSIONS'
1373 include 'COMMON.GEO'
1374 include 'COMMON.VAR'
1375 include 'COMMON.CHAIN'
1376 include 'COMMON.IOUNITS'
1377 include 'COMMON.CONTROL'
1378 include 'COMMON.LOCAL'
1379 include 'COMMON.INTERACT'
1380 c Read coordinates from input
1382 read(kanal,'(8f10.5)',end=10,err=10)
1383 & ((c(l,k),l=1,3),k=1,nres),
1384 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1387 c(j,2*nres)=c(j,nres)
1389 call int_from_cart1(.false.)
1392 dc(j,i)=c(j,i+1)-c(j,i)
1393 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1397 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1399 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1400 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1408 c----------------------------------------------------------------------------
1409 subroutine read_threadbase
1410 implicit real*8 (a-h,o-z)
1411 include 'DIMENSIONS'
1412 include 'COMMON.IOUNITS'
1413 include 'COMMON.GEO'
1414 include 'COMMON.VAR'
1415 include 'COMMON.INTERACT'
1416 include 'COMMON.LOCAL'
1417 include 'COMMON.NAMES'
1418 include 'COMMON.CHAIN'
1419 include 'COMMON.FFIELD'
1420 include 'COMMON.SBRIDGE'
1421 include 'COMMON.HEADER'
1422 include 'COMMON.CONTROL'
1423 include 'COMMON.DBASE'
1424 include 'COMMON.THREAD'
1425 include 'COMMON.TIME1'
1426 C Read pattern database for threading.
1427 read (icbase,*) nseq
1429 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1430 & nres_base(2,i),nres_base(3,i)
1431 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1433 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1434 c & nres_base(2,i),nres_base(3,i)
1435 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1439 if (weidis.eq.0.0D0) weidis=0.1D0
1448 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1449 write (iout,'(a,i5)') 'nexcl: ',nexcl
1450 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1453 c------------------------------------------------------------------------------
1454 subroutine setup_var
1455 implicit real*8 (a-h,o-z)
1456 include 'DIMENSIONS'
1457 include 'COMMON.IOUNITS'
1458 include 'COMMON.GEO'
1459 include 'COMMON.VAR'
1460 include 'COMMON.INTERACT'
1461 include 'COMMON.LOCAL'
1462 include 'COMMON.NAMES'
1463 include 'COMMON.CHAIN'
1464 include 'COMMON.FFIELD'
1465 include 'COMMON.SBRIDGE'
1466 include 'COMMON.HEADER'
1467 include 'COMMON.CONTROL'
1468 include 'COMMON.DBASE'
1469 include 'COMMON.THREAD'
1470 include 'COMMON.TIME1'
1471 C Set up variable list.
1477 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1479 ialph(i,1)=nvar+nside
1483 if (indphi.gt.0) then
1485 else if (indback.gt.0) then
1490 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1493 c----------------------------------------------------------------------------
1494 subroutine gen_dist_constr
1495 C Generate CA distance constraints.
1496 implicit real*8 (a-h,o-z)
1497 include 'DIMENSIONS'
1498 include 'COMMON.IOUNITS'
1499 include 'COMMON.GEO'
1500 include 'COMMON.VAR'
1501 include 'COMMON.INTERACT'
1502 include 'COMMON.LOCAL'
1503 include 'COMMON.NAMES'
1504 include 'COMMON.CHAIN'
1505 include 'COMMON.FFIELD'
1506 include 'COMMON.SBRIDGE'
1507 include 'COMMON.HEADER'
1508 include 'COMMON.CONTROL'
1509 include 'COMMON.DBASE'
1510 include 'COMMON.THREAD'
1511 include 'COMMON.TIME1'
1512 dimension itype_pdb(maxres)
1513 common /pizda/ itype_pdb
1515 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1516 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1517 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1519 do i=nstart_sup,nstart_sup+nsup-1
1520 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1521 cd & ' seq_pdb', restyp(itype_pdb(i))
1522 do j=i+2,nstart_sup+nsup-1
1524 ihpb(nhpb)=i+nstart_seq-nstart_sup
1525 jhpb(nhpb)=j+nstart_seq-nstart_sup
1527 dhpb(nhpb)=dist(i,j)
1530 cd write (iout,'(a)') 'Distance constraints:'
1535 cd if (ii.gt.nres) then
1540 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1541 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1542 cd & dhpb(i),forcon(i)
1546 c----------------------------------------------------------------------------
1548 implicit real*8 (a-h,o-z)
1549 include 'DIMENSIONS'
1550 include 'COMMON.MAP'
1551 include 'COMMON.IOUNITS'
1552 character*3 angid(4) /'THE','PHI','ALP','OME'/
1553 character*80 mapcard,ucase
1555 read (inp,'(a)') mapcard
1556 mapcard=ucase(mapcard)
1557 if (index(mapcard,'PHI').gt.0) then
1559 else if (index(mapcard,'THE').gt.0) then
1561 else if (index(mapcard,'ALP').gt.0) then
1563 else if (index(mapcard,'OME').gt.0) then
1566 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1567 stop 'Error - illegal variable spec in MAP card.'
1569 call readi (mapcard,'RES1',res1(imap),0)
1570 call readi (mapcard,'RES2',res2(imap),0)
1571 if (res1(imap).eq.0) then
1572 res1(imap)=res2(imap)
1573 else if (res2(imap).eq.0) then
1574 res2(imap)=res1(imap)
1576 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1578 & 'Error - illegal definition of variable group in MAP.'
1579 stop 'Error - illegal definition of variable group in MAP.'
1581 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1582 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1583 call readi(mapcard,'NSTEP',nstep(imap),0)
1584 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1586 & 'Illegal boundary and/or step size specification in MAP.'
1587 stop 'Illegal boundary and/or step size specification in MAP.'
1592 c----------------------------------------------------------------------------
1594 implicit real*8 (a-h,o-z)
1595 include 'DIMENSIONS'
1596 include 'COMMON.IOUNITS'
1597 include 'COMMON.GEO'
1598 include 'COMMON.CSA'
1599 include 'COMMON.BANK'
1600 include 'COMMON.CONTROL'
1602 character*620 mcmcard
1603 call card_concat(mcmcard)
1605 call readi(mcmcard,'NCONF',nconf,50)
1606 call readi(mcmcard,'NADD',nadd,0)
1607 call readi(mcmcard,'JSTART',jstart,1)
1608 call readi(mcmcard,'JEND',jend,1)
1609 call readi(mcmcard,'NSTMAX',nstmax,500000)
1610 call readi(mcmcard,'N0',n0,1)
1611 call readi(mcmcard,'N1',n1,6)
1612 call readi(mcmcard,'N2',n2,4)
1613 call readi(mcmcard,'N3',n3,0)
1614 call readi(mcmcard,'N4',n4,0)
1615 call readi(mcmcard,'N5',n5,0)
1616 call readi(mcmcard,'N6',n6,10)
1617 call readi(mcmcard,'N7',n7,0)
1618 call readi(mcmcard,'N8',n8,0)
1619 call readi(mcmcard,'N9',n9,0)
1620 call readi(mcmcard,'N14',n14,0)
1621 call readi(mcmcard,'N15',n15,0)
1622 call readi(mcmcard,'N16',n16,0)
1623 call readi(mcmcard,'N17',n17,0)
1624 call readi(mcmcard,'N18',n18,0)
1626 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1628 call readi(mcmcard,'NDIFF',ndiff,2)
1629 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1630 call readi(mcmcard,'IS1',is1,1)
1631 call readi(mcmcard,'IS2',is2,8)
1632 call readi(mcmcard,'NRAN0',nran0,4)
1633 call readi(mcmcard,'NRAN1',nran1,2)
1634 call readi(mcmcard,'IRR',irr,1)
1635 call readi(mcmcard,'NSEED',nseed,20)
1636 call readi(mcmcard,'NTOTAL',ntotal,10000)
1637 call reada(mcmcard,'CUT1',cut1,2.0d0)
1638 call reada(mcmcard,'CUT2',cut2,5.0d0)
1639 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1640 call readi(mcmcard,'ICMAX',icmax,3)
1641 call readi(mcmcard,'IRESTART',irestart,0)
1642 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1645 call reada(mcmcard,'DELE',dele,20.0d0)
1646 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1647 call readi(mcmcard,'IREF',iref,0)
1648 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1649 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1650 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1651 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1652 write (iout,*) "NCONF_IN",nconf_in
1655 c----------------------------------------------------------------------------
1656 cfmc subroutine mcmfread
1657 cfmc implicit real*8 (a-h,o-z)
1658 cfmc include 'DIMENSIONS'
1659 cfmc include 'COMMON.MCMF'
1660 cfmc include 'COMMON.IOUNITS'
1661 cfmc include 'COMMON.GEO'
1662 cfmc character*80 ucase
1663 cfmc character*620 mcmcard
1664 cfmc call card_concat(mcmcard)
1666 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1667 cfmc write(iout,*)'MAXRANT=',maxrant
1668 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1669 cfmc write(iout,*)'MAXFAM=',maxfam
1670 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1671 cfmc write(iout,*)'NNET1=',nnet1
1672 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1673 cfmc write(iout,*)'NNET2=',nnet2
1674 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1675 cfmc write(iout,*)'NNET3=',nnet3
1676 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1677 cfmc write(iout,*)'ILASTT=',ilastt
1678 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1679 cfmc write(iout,*)'MAXSTR=',maxstr
1680 cfmc maxstr_f=maxstr/maxfam
1681 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1682 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1683 cfmc write(iout,*)'NMCMF=',nmcmf
1684 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1685 cfmc write(iout,*)'IFOCUS=',ifocus
1686 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1687 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1688 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1689 cfmc write(iout,*)'INTPRT=',intprt
1690 cfmc call readi(mcmcard,'IPRT',iprt,100)
1691 cfmc write(iout,*)'IPRT=',iprt
1692 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1693 cfmc write(iout,*)'IMAXTR=',imaxtr
1694 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1695 cfmc write(iout,*)'MAXEVEN=',maxeven
1696 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1697 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1698 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1699 cfmc write(iout,*)'INIMIN=',inimin
1700 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1701 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1702 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1703 cfmc write(iout,*)'NTHREAD=',nthread
1704 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1705 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1706 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1707 cfmc write(iout,*)'MAXPERT=',maxpert
1708 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1709 cfmc write(iout,*)'IRMSD=',irmsd
1710 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1711 cfmc write(iout,*)'DENEMIN=',denemin
1712 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1713 cfmc write(iout,*)'RCUT1S=',rcut1s
1714 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1715 cfmc write(iout,*)'RCUT1E=',rcut1e
1716 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1717 cfmc write(iout,*)'RCUT2S=',rcut2s
1718 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1719 cfmc write(iout,*)'RCUT2E=',rcut2e
1720 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1721 cfmc write(iout,*)'DPERT1=',d_pert1
1722 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1723 cfmc write(iout,*)'DPERT1A=',d_pert1a
1724 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1725 cfmc write(iout,*)'DPERT2=',d_pert2
1726 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1727 cfmc write(iout,*)'DPERT2A=',d_pert2a
1728 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1729 cfmc write(iout,*)'DPERT2B=',d_pert2b
1730 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1731 cfmc write(iout,*)'DPERT2C=',d_pert2c
1732 cfmc d_pert1=deg2rad*d_pert1
1733 cfmc d_pert1a=deg2rad*d_pert1a
1734 cfmc d_pert2=deg2rad*d_pert2
1735 cfmc d_pert2a=deg2rad*d_pert2a
1736 cfmc d_pert2b=deg2rad*d_pert2b
1737 cfmc d_pert2c=deg2rad*d_pert2c
1738 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1739 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1740 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1741 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1742 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1743 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1744 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1745 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1746 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1747 cfmc write(iout,*)'RCUTINI=',rcutini
1748 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1749 cfmc write(iout,*)'GRAT=',grat
1750 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1751 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1755 c----------------------------------------------------------------------------
1757 implicit real*8 (a-h,o-z)
1758 include 'DIMENSIONS'
1759 include 'COMMON.MCM'
1760 include 'COMMON.MCE'
1761 include 'COMMON.IOUNITS'
1763 character*320 mcmcard
1764 call card_concat(mcmcard)
1765 call readi(mcmcard,'MAXACC',maxacc,100)
1766 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1767 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1768 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1769 call readi(mcmcard,'MAXREPM',maxrepm,200)
1770 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1771 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1772 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1773 call reada(mcmcard,'E_UP',e_up,5.0D0)
1774 call reada(mcmcard,'DELTE',delte,0.1D0)
1775 call readi(mcmcard,'NSWEEP',nsweep,5)
1776 call readi(mcmcard,'NSTEPH',nsteph,0)
1777 call readi(mcmcard,'NSTEPC',nstepc,0)
1778 call reada(mcmcard,'TMIN',tmin,298.0D0)
1779 call reada(mcmcard,'TMAX',tmax,298.0D0)
1780 call readi(mcmcard,'NWINDOW',nwindow,0)
1781 call readi(mcmcard,'PRINT_MC',print_mc,0)
1782 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1783 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1784 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1785 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1786 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1787 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1788 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1789 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1790 if (nwindow.gt.0) then
1791 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1793 winlen(i)=winend(i)-winstart(i)+1
1796 if (tmax.lt.tmin) tmax=tmin
1797 if (tmax.eq.tmin) then
1801 if (nstepc.gt.0 .and. nsteph.gt.0) then
1802 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1803 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1805 C Probabilities of different move types
1806 sumpro_type(0)=0.0D0
1807 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1808 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1809 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1810 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1811 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1812 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1813 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1815 print *,'i',i,' sumprotype',sumpro_type(i)
1816 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1817 print *,'i',i,' sumprotype',sumpro_type(i)
1821 c----------------------------------------------------------------------------
1822 subroutine read_minim
1823 implicit real*8 (a-h,o-z)
1824 include 'DIMENSIONS'
1825 include 'COMMON.MINIM'
1826 include 'COMMON.IOUNITS'
1828 character*320 minimcard
1829 call card_concat(minimcard)
1830 call readi(minimcard,'MAXMIN',maxmin,2000)
1831 call readi(minimcard,'MAXFUN',maxfun,5000)
1832 call readi(minimcard,'MINMIN',minmin,maxmin)
1833 call readi(minimcard,'MINFUN',minfun,maxmin)
1834 call reada(minimcard,'TOLF',tolf,1.0D-2)
1835 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1836 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1837 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1838 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1839 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1840 & 'Options in energy minimization:'
1841 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1842 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1843 & 'MinMin:',MinMin,' MinFun:',MinFun,
1844 & ' TolF:',TolF,' RTolF:',RTolF
1847 c----------------------------------------------------------------------------
1848 subroutine read_angles(kanal,*)
1849 implicit real*8 (a-h,o-z)
1850 include 'DIMENSIONS'
1851 include 'COMMON.GEO'
1852 include 'COMMON.VAR'
1853 include 'COMMON.CHAIN'
1854 include 'COMMON.IOUNITS'
1855 include 'COMMON.CONTROL'
1856 c Read angles from input
1858 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1859 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1860 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1861 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1864 c 9/7/01 avoid 180 deg valence angle
1865 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1867 theta(i)=deg2rad*theta(i)
1868 phi(i)=deg2rad*phi(i)
1869 alph(i)=deg2rad*alph(i)
1870 omeg(i)=deg2rad*omeg(i)
1875 c----------------------------------------------------------------------------
1876 subroutine reada(rekord,lancuch,wartosc,default)
1878 character*(*) rekord,lancuch
1879 double precision wartosc,default
1882 iread=index(rekord,lancuch)
1883 if (iread.eq.0) then
1887 iread=iread+ilen(lancuch)+1
1888 read (rekord(iread:),*,err=10,end=10) wartosc
1893 c----------------------------------------------------------------------------
1894 subroutine readi(rekord,lancuch,wartosc,default)
1896 character*(*) rekord,lancuch
1897 integer wartosc,default
1900 iread=index(rekord,lancuch)
1901 if (iread.eq.0) then
1905 iread=iread+ilen(lancuch)+1
1906 read (rekord(iread:),*,err=10,end=10) wartosc
1911 c----------------------------------------------------------------------------
1912 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1915 integer tablica(dim),default
1916 character*(*) rekord,lancuch
1923 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1924 if (iread.eq.0) return
1925 iread=iread+ilen(lancuch)+1
1926 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1929 c----------------------------------------------------------------------------
1930 subroutine multreada(rekord,lancuch,tablica,dim,default)
1933 double precision tablica(dim),default
1934 character*(*) rekord,lancuch
1941 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1942 if (iread.eq.0) return
1943 iread=iread+ilen(lancuch)+1
1944 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1947 c----------------------------------------------------------------------------
1948 subroutine openunits
1949 implicit real*8 (a-h,o-z)
1950 include 'DIMENSIONS'
1953 character*16 form,nodename
1956 include 'COMMON.SETUP'
1957 include 'COMMON.IOUNITS'
1959 include 'COMMON.CONTROL'
1960 integer lenpre,lenpot,ilen,lentmp
1962 character*3 out1file_text,ucase
1965 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1966 call getenv_loc("PREFIX",prefix)
1968 call getenv_loc("POT",pot)
1969 call getenv_loc("DIRTMP",tmpdir)
1970 call getenv_loc("CURDIR",curdir)
1971 call getenv_loc("OUT1FILE",out1file_text)
1972 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1973 out1file_text=ucase(out1file_text)
1974 if (out1file_text(1:1).eq."Y") then
1977 out1file=fg_rank.gt.0
1982 if (lentmp.gt.0) then
1983 write (*,'(80(1h!))')
1984 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1985 write (*,'(80(1h!))')
1986 write (*,*)"All output files will be on node /tmp directory."
1988 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1989 if (me.eq.king) then
1990 write (*,*) "The master node is ",nodename
1991 else if (fg_rank.eq.0) then
1992 write (*,*) "I am the CG slave node ",nodename
1994 write (*,*) "I am the FG slave node ",nodename
1997 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1998 lenpre = lentmp+lenpre+1
2000 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2001 C Get the names and open the input files
2002 #if defined(WINIFL) || defined(WINPGI)
2003 open(1,file=pref_orig(:ilen(pref_orig))//
2004 & '.inp',status='old',readonly,shared)
2005 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2006 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2007 C Get parameter filenames and open the parameter files.
2008 call getenv_loc('BONDPAR',bondname)
2009 open (ibond,file=bondname,status='old',readonly,shared)
2010 call getenv_loc('THETPAR',thetname)
2011 open (ithep,file=thetname,status='old',readonly,shared)
2012 call getenv_loc('ROTPAR',rotname)
2013 open (irotam,file=rotname,status='old',readonly,shared)
2014 call getenv_loc('TORPAR',torname)
2015 open (itorp,file=torname,status='old',readonly,shared)
2016 call getenv_loc('TORDPAR',tordname)
2017 open (itordp,file=tordname,status='old',readonly,shared)
2018 call getenv_loc('TORKCC',torkccname)
2019 open (itorkcc,file=torkccname,status='old',readonly,shared)
2020 call getenv_loc('THETKCC',thetkccname)
2021 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2022 call getenv_loc('FOURIER',fouriername)
2023 open (ifourier,file=fouriername,status='old',readonly,shared)
2024 call getenv_loc('ELEPAR',elename)
2025 open (ielep,file=elename,status='old',readonly,shared)
2026 call getenv_loc('SIDEPAR',sidename)
2027 open (isidep,file=sidename,status='old',readonly,shared)
2028 call getenv_loc('LIPTRANPAR',liptranname)
2029 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2030 #elif (defined CRAY) || (defined AIX)
2031 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2033 c print *,"Processor",myrank," opened file 1"
2034 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2035 c print *,"Processor",myrank," opened file 9"
2036 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2037 C Get parameter filenames and open the parameter files.
2038 call getenv_loc('BONDPAR',bondname)
2039 open (ibond,file=bondname,status='old',action='read')
2040 c print *,"Processor",myrank," opened file IBOND"
2041 call getenv_loc('THETPAR',thetname)
2042 open (ithep,file=thetname,status='old',action='read')
2043 c print *,"Processor",myrank," opened file ITHEP"
2044 call getenv_loc('ROTPAR',rotname)
2045 open (irotam,file=rotname,status='old',action='read')
2046 c print *,"Processor",myrank," opened file IROTAM"
2047 call getenv_loc('TORPAR',torname)
2048 open (itorp,file=torname,status='old',action='read')
2049 c print *,"Processor",myrank," opened file ITORP"
2050 call getenv_loc('TORDPAR',tordname)
2051 open (itordp,file=tordname,status='old',action='read')
2052 call getenv_loc('TORKCC',torkccname)
2053 open (itorkcc,file=torkccname,status='old',action='read')
2054 call getenv_loc('THETKCC',thetkccname)
2055 open (ithetkcc,file=thetkccname,status='old',action='read')
2056 c print *,"Processor",myrank," opened file ITORDP"
2057 call getenv_loc('SCCORPAR',sccorname)
2058 open (isccor,file=sccorname,status='old',action='read')
2059 c print *,"Processor",myrank," opened file ISCCOR"
2060 call getenv_loc('FOURIER',fouriername)
2061 open (ifourier,file=fouriername,status='old',action='read')
2062 c print *,"Processor",myrank," opened file IFOURIER"
2063 call getenv_loc('ELEPAR',elename)
2064 open (ielep,file=elename,status='old',action='read')
2065 c print *,"Processor",myrank," opened file IELEP"
2066 call getenv_loc('SIDEPAR',sidename)
2067 open (isidep,file=sidename,status='old',action='read')
2068 call getenv_loc('LIPTRANPAR',liptranname)
2069 open (iliptranpar,file=liptranname,status='old',action='read')
2070 c print *,"Processor",myrank," opened file ISIDEP"
2071 c print *,"Processor",myrank," opened parameter files"
2073 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2074 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2075 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2076 C Get parameter filenames and open the parameter files.
2077 call getenv_loc('BONDPAR',bondname)
2078 open (ibond,file=bondname,status='old')
2079 call getenv_loc('THETPAR',thetname)
2080 open (ithep,file=thetname,status='old')
2081 call getenv_loc('ROTPAR',rotname)
2082 open (irotam,file=rotname,status='old')
2083 call getenv_loc('TORPAR',torname)
2084 open (itorp,file=torname,status='old')
2085 call getenv_loc('TORDPAR',tordname)
2086 open (itordp,file=tordname,status='old')
2087 call getenv_loc('TORKCC',torkccname)
2088 open (itorkcc,file=torkccname,status='old')
2089 call getenv_loc('THETKCC',thetkccname)
2090 open (ithetkcc,file=thetkccname,status='old')
2091 call getenv_loc('SCCORPAR',sccorname)
2092 open (isccor,file=sccorname,status='old')
2093 call getenv_loc('FOURIER',fouriername)
2094 open (ifourier,file=fouriername,status='old')
2095 call getenv_loc('ELEPAR',elename)
2096 open (ielep,file=elename,status='old')
2097 call getenv_loc('SIDEPAR',sidename)
2098 open (isidep,file=sidename,status='old')
2099 call getenv_loc('LIPTRANPAR',liptranname)
2100 open (iliptranpar,file=liptranname,status='old')
2102 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2104 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2105 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2106 C Get parameter filenames and open the parameter files.
2107 call getenv_loc('BONDPAR',bondname)
2108 open (ibond,file=bondname,status='old',readonly)
2109 call getenv_loc('THETPAR',thetname)
2110 open (ithep,file=thetname,status='old',readonly)
2111 call getenv_loc('ROTPAR',rotname)
2112 open (irotam,file=rotname,status='old',readonly)
2113 call getenv_loc('TORPAR',torname)
2114 open (itorp,file=torname,status='old',readonly)
2115 call getenv_loc('TORDPAR',tordname)
2116 open (itordp,file=tordname,status='old',readonly)
2117 call getenv_loc('TORKCC',torkccname)
2118 open (itorkcc,file=torkccname,status='old',readonly)
2119 call getenv_loc('THETKCC',thetkccname)
2120 open (ithetkcc,file=thetkccname,status='old',readonly)
2121 call getenv_loc('SCCORPAR',sccorname)
2122 open (isccor,file=sccorname,status='old',readonly)
2124 call getenv_loc('THETPARPDB',thetname_pdb)
2125 print *,"thetname_pdb ",thetname_pdb
2126 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2127 print *,ithep_pdb," opened"
2129 call getenv_loc('FOURIER',fouriername)
2130 open (ifourier,file=fouriername,status='old',readonly)
2131 call getenv_loc('ELEPAR',elename)
2132 open (ielep,file=elename,status='old',readonly)
2133 call getenv_loc('SIDEPAR',sidename)
2134 open (isidep,file=sidename,status='old',readonly)
2135 call getenv_loc('LIPTRANPAR',liptranname)
2136 open (iliptranpar,file=liptranname,status='old',action='read')
2138 call getenv_loc('ROTPARPDB',rotname_pdb)
2139 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2142 call getenv_loc('TUBEPAR',tubename)
2143 open (itube,file=tubename,status='old',readonly)
2146 C 8/9/01 In the newest version SCp interaction constants are read from a file
2147 C Use -DOLDSCP to use hard-coded constants instead.
2149 call getenv_loc('SCPPAR',scpname)
2150 #if defined(WINIFL) || defined(WINPGI)
2151 open (iscpp,file=scpname,status='old',readonly,shared)
2152 #elif (defined CRAY) || (defined AIX)
2153 open (iscpp,file=scpname,status='old',action='read')
2155 open (iscpp,file=scpname,status='old')
2157 open (iscpp,file=scpname,status='old',readonly)
2160 call getenv_loc('PATTERN',patname)
2161 #if defined(WINIFL) || defined(WINPGI)
2162 open (icbase,file=patname,status='old',readonly,shared)
2163 #elif (defined CRAY) || (defined AIX)
2164 open (icbase,file=patname,status='old',action='read')
2166 open (icbase,file=patname,status='old')
2168 open (icbase,file=patname,status='old',readonly)
2171 C Open output file only for CG processes
2172 c print *,"Processor",myrank," fg_rank",fg_rank
2173 if (fg_rank.eq.0) then
2175 if (nodes.eq.1) then
2178 npos = dlog10(dfloat(nodes-1))+1
2180 if (npos.lt.3) npos=3
2181 write (liczba,'(i1)') npos
2182 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2184 write (liczba,form) me
2185 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2186 & liczba(:ilen(liczba))
2187 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2189 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2191 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2192 & liczba(:ilen(liczba))//'.mol2'
2193 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2194 & liczba(:ilen(liczba))//'.stat'
2196 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2197 & //liczba(:ilen(liczba))//'.stat')
2198 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2201 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2202 & liczba(:ilen(liczba))//'.const'
2207 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2208 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2209 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2210 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2211 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2213 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2215 rest2name=prefix(:ilen(prefix))//'.rst'
2217 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2220 #if defined(AIX) || defined(PGI)
2221 if (me.eq.king .or. .not. out1file)
2222 & open(iout,file=outname,status='unknown')
2224 if (fg_rank.gt.0) then
2225 write (liczba,'(i3.3)') myrank/nfgtasks
2226 write (ll,'(bz,i3.3)') fg_rank
2227 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2232 open(igeom,file=intname,status='unknown',position='append')
2233 open(ipdb,file=pdbname,status='unknown')
2234 open(imol2,file=mol2name,status='unknown')
2235 open(istat,file=statname,status='unknown',position='append')
2237 c1out open(iout,file=outname,status='unknown')
2240 if (me.eq.king .or. .not.out1file)
2241 & open(iout,file=outname,status='unknown')
2243 if (fg_rank.gt.0) then
2244 write (liczba,'(i3.3)') myrank/nfgtasks
2245 write (ll,'(bz,i3.3)') fg_rank
2246 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2251 open(igeom,file=intname,status='unknown',access='append')
2252 open(ipdb,file=pdbname,status='unknown')
2253 open(imol2,file=mol2name,status='unknown')
2254 open(istat,file=statname,status='unknown',access='append')
2256 c1out open(iout,file=outname,status='unknown')
2259 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2260 csa_seed=prefix(:lenpre)//'.CSA.seed'
2261 csa_history=prefix(:lenpre)//'.CSA.history'
2262 csa_bank=prefix(:lenpre)//'.CSA.bank'
2263 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2264 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2265 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2266 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2267 csa_int=prefix(:lenpre)//'.int'
2268 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2269 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2270 csa_in=prefix(:lenpre)//'.CSA.in'
2271 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2274 write (iout,'(80(1h-))')
2275 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2276 write (iout,'(80(1h-))')
2277 write (iout,*) "Input file : ",
2278 & pref_orig(:ilen(pref_orig))//'.inp'
2279 write (iout,*) "Output file : ",
2280 & outname(:ilen(outname))
2282 write (iout,*) "Sidechain potential file : ",
2283 & sidename(:ilen(sidename))
2285 write (iout,*) "SCp potential file : ",
2286 & scpname(:ilen(scpname))
2288 write (iout,*) "Electrostatic potential file : ",
2289 & elename(:ilen(elename))
2290 write (iout,*) "Cumulant coefficient file : ",
2291 & fouriername(:ilen(fouriername))
2292 write (iout,*) "Torsional parameter file : ",
2293 & torname(:ilen(torname))
2294 write (iout,*) "Double torsional parameter file : ",
2295 & tordname(:ilen(tordname))
2296 write (iout,*) "SCCOR parameter file : ",
2297 & sccorname(:ilen(sccorname))
2298 write (iout,*) "Bond & inertia constant file : ",
2299 & bondname(:ilen(bondname))
2300 write (iout,*) "Bending parameter file : ",
2301 & thetname(:ilen(thetname))
2302 write (iout,*) "Rotamer parameter file : ",
2303 & rotname(:ilen(rotname))
2304 write (iout,*) "Thetpdb parameter file : ",
2305 & thetname_pdb(:ilen(thetname_pdb))
2306 write (iout,*) "Threading database : ",
2307 & patname(:ilen(patname))
2309 &write (iout,*)" DIRTMP : ",
2311 write (iout,'(80(1h-))')
2315 c----------------------------------------------------------------------------
2316 subroutine card_concat(card)
2317 implicit real*8 (a-h,o-z)
2318 include 'DIMENSIONS'
2319 include 'COMMON.IOUNITS'
2321 character*80 karta,ucase
2323 read (inp,'(a)') karta
2326 do while (karta(80:80).eq.'&')
2327 card=card(:ilen(card)+1)//karta(:79)
2328 read (inp,'(a)') karta
2331 card=card(:ilen(card)+1)//karta
2334 c----------------------------------------------------------------------------------
2336 implicit real*8 (a-h,o-z)
2337 include 'DIMENSIONS'
2338 include 'COMMON.CHAIN'
2339 include 'COMMON.IOUNITS'
2341 open(irest2,file=rest2name,status='unknown')
2342 read(irest2,*) totT,EK,potE,totE,t_bath
2345 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2348 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2351 read (irest2,*) iset
2356 c---------------------------------------------------------------------------------
2357 subroutine read_fragments
2358 implicit real*8 (a-h,o-z)
2359 include 'DIMENSIONS'
2363 include 'COMMON.SETUP'
2364 include 'COMMON.CHAIN'
2365 include 'COMMON.IOUNITS'
2367 include 'COMMON.CONTROL'
2368 read(inp,*) nset,nfrag,npair,nfrag_back
2369 if(me.eq.king.or..not.out1file)
2370 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2371 & " nfrag_back",nfrag_back
2373 read(inp,*) mset(iset)
2375 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2377 if(me.eq.king.or..not.out1file)
2378 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2379 & ifrag(2,i,iset), qinfrag(i,iset)
2382 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2384 if(me.eq.king.or..not.out1file)
2385 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2386 & ipair(2,i,iset), qinpair(i,iset)
2389 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2390 & wfrag_back(3,i,iset),
2391 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2392 if(me.eq.king.or..not.out1file)
2393 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2394 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2399 C---------------------------------------------------------------------------
2400 subroutine read_afminp
2401 implicit real*8 (a-h,o-z)
2402 include 'DIMENSIONS'
2406 include 'COMMON.SETUP'
2407 include 'COMMON.CONTROL'
2408 include 'COMMON.CHAIN'
2409 include 'COMMON.IOUNITS'
2410 include 'COMMON.SBRIDGE'
2411 character*320 afmcard
2413 call card_concat(afmcard)
2414 call readi(afmcard,"BEG",afmbeg,0)
2415 call readi(afmcard,"END",afmend,0)
2416 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2417 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2418 print *,'FORCE=' ,forceAFMconst
2419 CCCC NOW PROPERTIES FOR AFM
2422 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2424 distafminit=dsqrt(distafminit)
2425 print *,'initdist',distafminit
2429 c-------------------------------------------------------------------------------
2430 subroutine read_dist_constr
2431 implicit real*8 (a-h,o-z)
2432 include 'DIMENSIONS'
2436 include 'COMMON.SETUP'
2437 include 'COMMON.CONTROL'
2438 include 'COMMON.CHAIN'
2439 include 'COMMON.IOUNITS'
2440 include 'COMMON.SBRIDGE'
2441 integer ifrag_(2,100),ipair_(2,100)
2442 double precision wfrag_(100),wpair_(100)
2443 character*500 controlcard
2445 write (iout,*) "Calling read_dist_constr"
2446 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2448 call card_concat(controlcard)
2449 call readi(controlcard,"NFRAG",nfrag_,0)
2450 call readi(controlcard,"NPAIR",npair_,0)
2451 call readi(controlcard,"NDIST",ndist_,0)
2452 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2453 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2454 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2455 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2456 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2457 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2458 c write (iout,*) "IFRAG"
2460 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2462 c write (iout,*) "IPAIR"
2464 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2468 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2469 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2470 & ifrag_(2,i)=nstart_sup+nsup-1
2471 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2473 if (wfrag_(i).gt.0.0d0) then
2474 do j=ifrag_(1,i),ifrag_(2,i)-1
2475 do k=j+1,ifrag_(2,i)
2476 c write (iout,*) "j",j," k",k
2478 if (constr_dist.eq.1) then
2483 forcon(nhpb)=wfrag_(i)
2484 else if (constr_dist.eq.2) then
2485 if (ddjk.le.dist_cut) then
2490 forcon(nhpb)=wfrag_(i)
2497 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2500 if (.not.out1file .or. me.eq.king)
2501 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2502 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2504 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2505 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2512 if (wpair_(i).gt.0.0d0) then
2520 do j=ifrag_(1,ii),ifrag_(2,ii)
2521 do k=ifrag_(1,jj),ifrag_(2,jj)
2525 forcon(nhpb)=wpair_(i)
2526 dhpb(nhpb)=dist(j,k)
2528 if (.not.out1file .or. me.eq.king)
2529 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2530 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2532 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2533 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2541 if (constr_dist.eq.11) then
2542 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2543 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2544 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2547 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2548 & ibecarb(i),forcon(nhpb+1)
2550 if (forcon(nhpb+1).gt.0.0d0) then
2552 if (ibecarb(i).gt.0) then
2553 ihpb(i)=ihpb(i)+nres
2554 jhpb(i)=jhpb(i)+nres
2556 if (dhpb(nhpb).eq.0.0d0)
2557 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2559 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2560 C if (forcon(nhpb+1).gt.0.0d0) then
2562 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2564 if (.not.out1file .or. me.eq.king)
2565 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2566 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2568 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2569 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2576 c-------------------------------------------------------------------------------
2578 subroutine flush(iu)
2583 subroutine flush(iu)
2588 c------------------------------------------------------------------------------
2589 subroutine copy_to_tmp(source)
2590 include "DIMENSIONS"
2591 include "COMMON.IOUNITS"
2592 character*(*) source
2593 character* 256 tmpfile
2597 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2598 inquire(file=tmpfile,exist=ex)
2600 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2601 & " to temporary directory..."
2602 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2603 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2607 c------------------------------------------------------------------------------
2608 subroutine move_from_tmp(source)
2609 include "DIMENSIONS"
2610 include "COMMON.IOUNITS"
2611 character*(*) source
2614 write (*,*) "Moving ",source(:ilen(source)),
2615 & " from temporary directory to working directory"
2616 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2617 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2620 c------------------------------------------------------------------------------
2621 subroutine random_init(seed)
2623 C Initialize random number generator
2625 implicit real*8 (a-h,o-z)
2626 include 'DIMENSIONS'
2629 logical OKRandom, prng_restart
2631 integer iseed_array(4)
2633 include 'COMMON.IOUNITS'
2634 include 'COMMON.TIME1'
2635 include 'COMMON.THREAD'
2636 include 'COMMON.SBRIDGE'
2637 include 'COMMON.CONTROL'
2638 include 'COMMON.MCM'
2639 include 'COMMON.MAP'
2640 include 'COMMON.HEADER'
2641 include 'COMMON.CSA'
2642 include 'COMMON.CHAIN'
2643 include 'COMMON.MUCA'
2645 include 'COMMON.FFIELD'
2646 include 'COMMON.SETUP'
2647 iseed=-dint(dabs(seed))
2648 if (iseed.eq.0) then
2649 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2650 & 'Random seed undefined. The program will stop.'
2651 write (*,'(/80(1h*)/20x,a/80(1h*))')
2652 & 'Random seed undefined. The program will stop.'
2654 call mpi_finalize(mpi_comm_world,ierr)
2656 stop 'Bad random seed.'
2659 if (fg_rank.eq.0) then
2663 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2664 OKRandom = prng_restart(me,iseed)
2667 tmp=65536.0d0**(4-i)
2668 iseed_array(i) = dint(seed/tmp)
2669 seed=seed-iseed_array(i)*tmp
2672 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2673 & (iseed_array(i),i=1,4)
2674 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2675 & (iseed_array(i),i=1,4)
2676 OKRandom = prng_restart(me,iseed_array)
2679 c r1 = prng_next(me)
2680 r1=ran_number(0.0D0,1.0D0)
2682 & write (iout,*) 'ran_num',r1
2683 if (r1.lt.0.0d0) OKRandom=.false.
2685 if (.not.OKRandom) then
2686 write (iout,*) 'PRNG IS NOT WORKING!!!'
2687 print *,'PRNG IS NOT WORKING!!!'
2690 call mpi_abort(mpi_comm_world,error_msg,ierr)
2693 write (iout,*) 'too many processors for parallel prng'
2694 write (*,*) 'too many processors for parallel prng'
2702 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)